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ABSTRACT 

We develop and implement numerical methods for including stellar feedback in galaxy-scale numerical 
simulations. Our models include simplified treatments of heating by Type-I and Type-II supernovae, gas 
recycling from young stars and AGB winds, heating from the shocked stellar winds, HII photoionization 
heating, and radiation pressure from stellar photons. The energetics and time-dependence associated with 
the feedback are taken directly from stellar evolution models. We implement these stellar feedback models 
in smoothed particle hydrodynamic simulations with pc-scale resolution, modeling galaxies from SMC-like 
dwarfs and Milky-way analogues to massive z ^ 2 star forming disks. Absent stellar feedback, gas cools 
rapidly and collapses without limit into dense sub-units, inconsistent with observations. By contrast, in all 
cases with feedback, the interstellar medium (ISM) quickly approaches a statistical steady state in which 
giant molecular clouds (GMCs) continuously form, disperse, and re-form, leading to a multi-phase ISM. In 
this paper, we quantify the properties of the ISM and GMCs in this self -regulated state. In a companion paper 
we study the galactic winds driven by stellar feedback. 

Our primary results on the structure of the ISM in star forming galaxies include: 

1) Star forming galaxies generically self-regulate so that the cool, dense gas maintains Toomre's 2^1- 
Most of the volume is occupied by relatively diffuse hot gas, while most of the mass is in dense GMC 
complexes created by self-gravity. The phase structure of the gas and the gas mass fraction at high densities 
are much more sensitive probes of the physics of stellar feedback than integrated quantities such as the 
Toomre Q or gas velocity dispersion. 

2) Different stellar feedback mechanisms act on different spatial (and density) scales. Radiation pressure 
and HII gas pressure are critical for preventing runaway collapse of dense gas in GMCs. Shocked supernovae 
ejecta and stellar winds dominate the dynamics of the volume-filling hot gas. However, this gas primarily 
vents out of the star-forming disk and contributes only modestly to the mid-plane ISM pressure. 

3) The galaxy-averaged star formation rate is determined by feedback, with different mechanisms dom- 
inating in different galaxy types. For a given feedback efficiency, restricting star formation to molecular gas 
or modifying the cooling function has little effect on the star formation rate in the galaxies we model (includ- 
ing an SMC-mass dwarf). By contrast, changing the feedback mechanisms or assumed feedback efficiencies 
directly translates to shifts off of the observed Kennicutt-Schmidt relation. 

4) Self-gravity leads to GMCs with an approximately self-similar mass function cx M~^, with a high 
mass cutoff determined by the characteristic Jeans/Toomre mass of the system. In all of our galaxy models, 
GMCs live for a few dynamical times before they are disrupted by stellar feedback. The net star formation 
efficiency in GMCs ranges from ^ 1 % in dwarfs and MW-like spirals to nearly ^ 10% in gas-rich rapidly star 
forming galaxies. GMCs are approximately virialized, but there is a large dispersion in the virial parameter 
for a given GMC mass and lower mass GMCs tend to be preferentially unbound. 

Key words: galaxies: formation — star formation: general — galaxies: evolution — galaxies: active — 
cosmology: theory 
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1 INTRODUCTION 

The Kennicutt-Schmidt (KS) law for star formation in galaxies im- 
plies a gas consumption time of ~ 50 dynamical times ( [Kennicuttl 
|1998[ >; comparably long gas consumption times are inferred even 
on the scale of giant molecular clouds (GMCs) and other dense 
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Table 1. Galaxy Models 
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Parameters describing our galaxy models, used as the initial conditions for all of the simulations: 

(1) Model name: shorthand for models of a high-redshift massive starburst (HiZ); local gas-rich galaxy (Sbc); MW-analogue (MW); and isolated 
SMC-mass dwarf (SMC). (2) e^: gravitational force softening in our highest-resolution simulations. (3) m, : gas particle mass in our highest-resolution 
simulations. New star particles formed have mass = 0.5 m, , disk/bulge particles m,, and dark matter halo particles 5 m, . (4) Mhaio: halo mass. (5) 
c: halo concentration. Values lie on the halo mass-concentration relation at the appropriate redshift (z = for SMC, Sbc, and MW; z = 2 for HiZ). (6) 
Vmax- halo maximum circular velocity. (7) Mbary: total baryonic mass. (8) Mb: bulge mass. (9) a: |Hemquist|{T990) profile scale-length for bulge. (10) 
Mj: stellar disk mass. (11) rj: stellar disk scale length. (12) h: stellar disk scale-height. (13) Mg: gas disk mass. (14) fg: gas disk scale length (gas 
scale-height determined so that Q= \). (15) /ga,: average gas fraction of the disk inside of the stellar Rg (Mg[< Sf,]/(Mg[< ] +M[i[< ReX)). The 
total gas fraction, including the extended disk, is ~ 50% larger. (15) ;jy„: gas disk dynamical time at the half-gas mass radius. 



star-forming regions (e.g. |Krumholz & Tan||2007^ . IVloreover, the 
total fraction of the gas turned into stars in GMCs over their life- 
time is only a few percent (Zuckerman & Evans 1974 , Wi lliams &| 
|McKee|199"7{|Evans|1999(,Evans et al.,2009^ . By contrast, absent 
stellar feedback, GMCs experience runaway collapse to densities 
much higher than observed, and eventually turn a near-unity frac- 
tion of their gas i nto stars ([Hopkins e t al.|201 la||Tasker|201 1 [|Bour-| 
riaud et al. 201 ()| | Dobbs et al.|2011[ ]Krumholz et al. 2011, Harper- 
Clark & Murrayj2011) . An analogous problem occurs in cosmo- 



logical models of galaxy evolution without strong stellar feedback: 
gas cools efficiently, collapsing into clumps that runaway to high 
densities on a few dynamical times and produce an order of magni- 
tude more stars than is observed (e.g. |Katz et al.|1996[|Somerville| 
|& Primackl 19991 [Cole et al.|2000| and references therein). 

Neither thermal pressure nor turbulence in and of itself can 
stave off collapse in the ISM: cooling is rapid, so that thermal sup- 
port is quickly lost, and turbulent support dissipates on a single 
crossing time (e.g., [Ostriker et aL||2001| l. Some mechanism must 
therefore continuously inject energy and momentum into the gas 
on both GMC and galactic scales, in order to prevent runaway col- 
lapse to arbitrarily high densities. Moreover, it is not enough to 
simply slow down the collapse of gas; if that were the case, the gas 
would eventually turn into stars on long timescales. Instead, mass 
must actually be expelled from both GMCs and galaxies as a whole. 

A number of physical mechanisms have been proposed as 
sources of random motions in the ISM and GMCs: these include 
photo-ionization, stellar winds, radiation pressure from both UV 
and IR photons, proto-stellar jets, cosmic-rays created in supernova 
shocks, and supernovae themselves (this is not an exhaustive list!). 
An alternative class of models suggests that the turbulence is not 
explicitly produced by star formation, but comes from gravitational 
cascades from larger scales or instabilities that can tap into the dif- 
ferential rotation in galaxies, (see e.g. [Mac Low & K lessen 2 004| 
[McKee & Ostriker 2007 , and references therein, for reviews of ISM 
turbulence and stellar feedback). 

Making a realistic comparison between models of the ISM 
and observations ultimately requires numerical simulations, since 
the systems of interest are highly non-linear, chaotic, and time- 
dependent, with a large number of competing physical processes 
and a huge dynamic range in space and time. However, most nu- 
merical experiments have traditionally been restricted to ideal- 
ized simulations of single star-forming regions or sub-regions of 
single GMCs. There are good reasons for doing so, since these 



calculations are critical for studying the physics of star forma- 
tion. They cannot, however, address a large number of observa- 
tional constraints such as the balance of phases in the ISM; the 
origin of the global Schmidt-Kennicutt relation; and observations 
of GMC/galaxy scaling relations, such as mass functions and the 
linewidth-size relation. Likewise, the lifetimes of GMCs and their 
internal structure cannot be addressed in detail without a model that 
self-consistently incorporates continuous GMC formation, dissolu- 
tion, and re-formation. In short, addressing these questions requires 
galaxy-scale simulations that also self-consistently follow feedback 
from sub-GMC through super-galactic scales. 

There are at least two major challenges that numerical sim- 
ulations must confront to model feedback and GMC-level phase 
structure explicitly. First, simulations must actually resolve the for- 
mation of GMCs and the phase structure of the ISM (spatial res- 
olution ~ Ipc; although in massive gas-rich galaxies GMCs can 
be much larger, ~kpc in size). Numerical simulations of isolated 
galaxies, galaxy mergers, and cosmological "zoom-in" simulations 
of individual halos can now reach this resolution, but it remains 
prohibitive in full cosmological simulations. 

The second challenge facing simulations is that a wide range 
of processes appear to be important sources of feedback and turbu- 
lence in the ISM of galaxies: e.g., protostellar jets, HII regions, stel- 
lar winds, radiation pressure from young stars, and supernovae. It is 
therefore unlikely that the full problem, which spans many orders 
of magnitude in spatial, mass, and time scales, can be treated with 
just a single simple prescription, without including multiple physi- 
cal processes. In particular the physics of the dense ISM is clearly 
very different from that of the diffuse ISM, and different feedback 
processes likely dominate in each case. Thus far, however, many 
calculations that model stellar feedback explicitly have invoked 



only thermal heating via supernovae (e.g. IThacker & Couchman 
2000| [Govemato et al.|[2007) ICev erino & Kly pin 2009 Teyss^ 
20 10 , ^la-Reese et al. 2011 , Pilkington 



et al.|20I0l|Coli'n et al. 



et al.|201 1 1. The thermal energy, however, tends to be rapidly radi- 



ated away, especially in dense sub-regions of GMCs. Perhaps even 
more importantly, the actual time required for SNe to explode (a 
few Myr) is longer ^an the observed "lifetime" of low-mass GMCs 
and so supernovae cannot be the dominant mechanism unbinding 
GMCs in local galaxies jEvans et al.|2009| >. 

In [Hopkins et al.| ( |20Ila^ (hereafter Paper 1) we developed a 
new numerical model for the effects of the momentum injection 
from massive stars on small scales in star-forming regions. Us- 
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ing galaxy-scale simulations with pc-scale resolution to resolve the 
formation of individual GMCs, we accounted (approximately) for 
the momentum deposited as gas and dust in the vicinity of young 
stars absorbs the stellar radiation. We showed that this feedback 
mechanism alone can dramatically suppress runaway star forma- 
tion, disrupting collapsing GMCs once a few percent of their mass 
is turned into stars, and naturally producing galaxies on the ob- 
served Schmidt-Kennicutt relation. This is, however, just a first step 
towards a more complete model of the ISM. 

In this paper, we therefore extend the models in Paper I to 
include HII photo-ionization heating and the momentum injec- 
tion, gas heating (energy injection) and mass recycling from su- 
pernovae (Types I and II) and stellar winds (both early fast winds 
and late, AGB winds)Q We also extend the radiation pressure 
model to include not just the small-scale radiation pressure force 
from multiply-scattered IR photons in GMCs, but also the long- 
range radiation pressure force produced by the diffuse interstellar 
UV through IR radiation field. Our models of these feedback pro- 
cesses are physically motivated and tied to stellar evolution calcu- 
lations; they are, however, also necessarily simplified and approx- 
imate given the range of spatial scales and processes we consider. 
In addition to an improved treatment of stellar feedback, we also 
extend our cooling and star formation models to include explicit 
models for molecular chemistry, allowing for the suppression of 
star formation via photo-dissociation and inefficient molecule for- 
mation. These are all processes that have been discussed analyti- 
cally in the literature, and in some cases explored in idealized sim- 
ulations that focus on an individual physical effect. However, by 
including all of these processes together, our study allows us to 
assess the relative importance of different feedback mechanisms 
and the interplay between the different feedback and star forma- 
tion processes in a given galaxy model. Towards this end, we have 
implemented this star formation and feedback physics in a range 
of idealized, observationally-motivated galaxy models that span 
low-redshift dwarfs. Milky Way-like systems, and nearby starbursts 
through to extremely massive, high-redshift starburst galaxies. 

The remainder of this paper is organized as follows. In ^we 
describe the galaxy models and our method of implementing feed- 
back from the various mechanisms outlined above. The Appendix 
contains tests of the numerical methods. In § |3] we discuss how 
each feedback mechanism singly and in combination effects the 
observed galaxy morphologies in both gas and stars, for a range of 
different galaxy properties. In § [4] we consider the star formation 
histories of each model and show how they are influenced by each 
feedback mechanism separately and in concert. We investigate the 
role of molecular chemistry in §[5] and show that it has little effect 
on galaxy properties or global star formation (for the galaxy models 
we consider). We examine the global properties of the ISM in §|6] 
including gas disk velocity dispersions, scale heights, and Toomre 
parameters, and we illustrate how different feedback mechanisms 
contribute to supporting turbulence and maintaining marginal disk 
stability. We consider the phase structure of the ISM in § [7] and 
show how the presence of these feedback mechanisms together 
naturally gives rise to a multi-phase ISM (with cold molecular, 
warm ionized, and hot diffuse components). In § [S] we examine 
the properties of the cold dense gas in the simulated GMC pop- 
ulation: GMC mass functions, densities, velocity dispersions and 
scaling laws, virial parameters, lifetimes, and star formation effi- 




' Movies of these simulations are available at 'h ttps : / / www . cf a ■ | 

|harvard . edu/-phopkins /Site/Research . html| 



Figure 1. Morphology of the gas and stars in a standard simulation of a HiZ 
disk model, a massive, z ~ 2 — 4 starburst disk with M* > 1 00 Mq yr~', 
with all feedback mechanisms enabled. The images are shown at ~ 2 or- 
bital times after the beginning of the simulation when the disk is in a 
feedback-regulated steady-state. We show face-on (left) and edge-on right 
projections. Top: Stars. The image is a mock ugr (SDSS-band) compos- 
ite, with the spectrum of all stars calculated from their known age and 
metallicity, and dust extinction/reddening accounted for from the line-of- 
sight dust mass. The brightness follows a logarithmic scale with a stretch of 
ft! 2dex. Bottom: Gas. Brightness encodes projected gas density (logarith- 
mically scaled with a ft 4dex stretch); color encodes gas temperature with 
blue being T < lOOOK molecular gas, pink ~ lO'* - 10^ K warm ionized 
gas, and yellow > 10* K hot gas. Gravitational collapse forms massive kpc- 
scale star cluster complexes that give rise to the clumpy morphology (edge 
on, similar to "clump-chain" galaxies) that is further enhanced in the optical 
by patchy extinction. Violent outflows are present, emerging both from the 
complexes and the disk as a whole, driven by the massive stai'burst. These 
outflows maintain a thick gas disk and dismpt many of the clumps, which 
continuously re-form. 



ciencies. We compare each of these properties to observations and 
show how they are affected by the presence or absence of different 
feedback mechanisms. Finally, in ^we summarize our results and 
discuss their implications. 

These results presented in this paper provide an explanation 
for disruption of GMCs and the observed slow rate of star for- 
mation in galaxies, but they do not address the question of why 
galaxies have a baryon fraction well below that of the universe as 
a whole. In a companion paper we show that the combined effects 
of radiation pressure and supemovae drive powerful galaxy scale 
winds, ejecting baryons from star forming galaxies. 

2 THE SIMULATIONS 

The initial conditions and numerical details of the simulations used 
here are described in detail in Paper I (see their Section 2 and Ta- 
bles 1-3). However we briefly summarize some basic properties of 
the models here. The simulations were performed with the paral- 
lel Tree-SPH code GADGET-3 jSpringel|2005l l. They include stars, 
dark matter, and gas, with cooling, shocks, star formation, and stel- 
lar feedback. 
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Table 2. Feedback Physics for Simulations Used in This Paper 
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Summary of the physics included or excluded in the simulations shown throughout this paper The implementation of each set of 
physics is described in §|2] A / means the given "module" was included, X means it was excluded. Mod="Modified" runs are 
described in the text. The Enforced SED models modify the calculation of the long-range radiation pressure and are discussed in the 
Appendix. Paper I simulations use the SF model of Molecular Model A. 




Figure 2. Stellar and gas images as in Figure[T] but for the She model, a 
dwarf starburst galaxy. The disk is again clumpy, but the clumps are rela- 
tively much smaller than those in the HiZ model, and do not gravitationally 
dominate - the global disk structure more clearly traces spiral structure al- 
beit with a disturbed, filamentary and irregular pattern in the gas. There is a 
stellai' bar in the central few kpc and a bright starburst, but the disk is quite 
dusty, especially edge-on. Strong outflows arise from the central few kpc, 
producing a clumpy, multi-phase wind. 



Gas follows a standard atomic cooling curve but in addition 
can cool to < 100 K via fine-structure cooling. This allows it to 
collapse to very high densities, and star formation occurs in dense 
regions above a threshold n > lOOOcm^'', with a rate p* — epjts 
where is the free-fall time and e = 1.5% is an efficiency taken 
from observations of star-forming regions with the same densities 
( [Knimholz & Tan|2007[ and references therein). In Paper 1 we show 




Figure 3. Stellar and gas images as in Figure [T] but for the Milky Way 
model. These images are for a model with a slightly higher dark matter 
fraction than our standard MW model, such that a bar does not form. The 
resulting morphology in gas and stars is a canonical spiral (compare e.g. 
MlOl), with dust lanes along the spiral arms, flocculant morphology and 
features such as spurs/feathers in the outer arms. The star formation is still 
concentrated in star clusters (visible as small blue knots, largely along the 
arms) - but these are relatively much smaller than the giant clumps in Fig- 
ure[T] owing to the much smaller Jeans mass. SNe-driven bubbles/holes are 
clear in the outer gas disk. 

that the equilibrium structure and SFR are basically independent 
of the small-scale star formation law (robust to order-of-magnitude 
variations in the threshold and efficiency, and changes in the power- 
law index p cc p'^") because it is feedback-regulated (so clumps 
simply form stars until sufficient feedback halts further collapse). 
The feedback models are implemented in four distinct initial 
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3 kpc 



Figure 4. Stellar and gas images for our standard Milky Way model after 
it develops a strong bar (compare with Fig. [3] which is our MW model for 
a somewhat higher dark matter halo mass, which does not develop a bar). 
The morphology is that of a typical barred spiral (compare e.g. NGC1097). 
Star cluster formation is enhanced along the bar, with a small nuclear con- 
centration of star formation driven by inflows along the bar. 




Figure 5. Stellar and gas images as in Figure [T] but for the SMC model, 
an isolated SMC-mass dwarf galaxy. The moiphology is here that of an 
in'egular galaxy, with a puffy/thick disk of mostly ionized gas, irregularly 
distributed star formation, and large SNe bubbles driven by overlapping 
SNe explosions. Despite a very low absolute star formation rate, a wind 
is plainly evident with a mix of hot gas and entrained cold/warm material, 
with various filaments/loops/arcs driven by the hot gas bubbles as well as 
directly venting hot gas. 



disk models, designed to represent a range of characteristic galaxy 
types. Each initial disk has a bulge, stellar and gaseous disk, halo, 
and central BH (although to isolate the role of stellar feedback, 
models for BH growth and feedback are disabled); they are initial- 
ized in equilibrium so that in the absence of cooling, star formation, 
and feedback there are no significant transients. The gaseous disk is 
initially vertically pressure-supported, but this cools away in much 
less than a dynamical time and the emergent vertical structure de- 
pends on feedback (changing the initial support to turbulent disper- 
sion does not change our results). Our "low" resolution runs (used 
to evolve the simulations for several Gyr, to ensure steady-state be- 
havior) use ~ 3 X 10* particles, with ~ 10* particles in the disk, 
giving SPH smoothing lengths of ~ lOpc in the central few kpc of 
a MW-like disk (scaling linearly with the disk size/mass scale). Our 
"standard" resolution cases use ~ 30 times as many particles, and 
correspondingly reach up to ~ 1 — 5 pc smoothing lengths and par- 
ticle masses of 500Mq, and are run for a few orbital times each. 
A couple of ultra-high resolution runs used for convergence tests 
employ ~ lO' particles, with sub-pc resolution on kpc scales. 
The disk models include: 

(1) SMC: an SMC-like dwarf, with baryonic mass Mba, = 
8.9 X \& M(V) and halo ma ss Mhaio = 2 x lO'" Mq (concentration 
c = 15), a Hernquist| i 1990 1 profile bulge with mass = 10^ Mq, 
and exponential stellar {mj = 1.3 x 10*Mq) and gas disks (mg — 
7.5 X 10** Mq) with scale-lengths hd = 0.7 and hg = 2.1 kpc, re- 
spectively. Table 1 gives the maximum circular velocity. The initial 
scale-height is zo ~ 140 pc and the stellar disk is initialized such 
that the Toomre Q — I everywhere. 

(2) MW: a MW-like galaxy, with halo (Mhaio, f) = 
(1.6 X IO'^Mq, 12), and baryonic (Mbar, mt, nid, nig) = 
(7.1,1.5,4.7,0.9) X 10^° Mq with scale-lengths {hd,hg,zo) = 
(3.0,6.0, 0.3) kpc. 

(3) Sbc: a LIRG-like galaxy (i.e. a more bary on-dominated 
gas-rich spiral characteristic of those observed at low red- 
shifts), with halo (Mhaio, c) = (1.5 x 10"Mq, 11), and baryonic 
{Mbm-,nib,md,mg) = (10.5,1.0,4.0,5.5) x IO'Mq with scale- 
lengths {hd, hg, zo) ~ (1-3, 2.6, 0.13)kpc. 

(4) HiZ: a high-redshift massive starburst disk, chosen to 
match the properties of the observed non-merging but mas- 
sively star-forming SMG population, with halo (Mhaio, c) = (1.4 x 
IO'^Mq, 3.5) and a virial radius appropriately rescaled for a halo 
at z = 2 rather than z = 0, and baryonic (Mbar, nii,, nid, irig) — 
(10.7,0.7,3,7) X IO'^Mq with scale-lengths {hd,hg,zo) = 
(1.6, 3.2,0.32)kpc. 

The salient galaxy properties are summarized in Table [T] re- 
produced from Paper I. The feedback and ISM physics we include 
in different runs is summarized in Table[2l 



2.1 Stellar Feedback 

2.1.1 Local Momentum Flux from Radiation Pressure, 
Supernovae, & Stellar Winds 

First, we consider the radiation pressure produced on a sub-GMC 
scale from young star clusters. This feedback mechanism is that 
presented and discussed extensively in Paper I, but we briefly sum- 
marize its properties here. At each timestep, gas particles identify 
the nearest density peak (via a friends-of-friends algorithm aver- 
aged over a scale of a couple smoothing lengths), representing the 
center of the nearest star-forming "clump" or GMC-analog. The 
gas then lies at some point Rd from the cloud center. The summed 
luminosity from the simulation star particles inside the radius L(< 
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Rci) is determined from each star's Imown age and metallicity ac- 
cording to tiie STARBURST99 model ( |Leitherer et al.|1999^ . As- 
suming tliis is equally distributed to all gas inside the cloud, the 
rate of momentum deposition is then just Piad ~ (1 +TiR)Lincident/c, 
where the term 1 + tir accounts for the fact that most of the ini- 
tial optical/UV radiation is absorbed and re-radiated in the IR; 
riR = Egas «;iR is the optical depth in the IR, which allows for the 
fact that the momentum is boosted by multiple scatterings in op- 
tically thick regions. Here Egas is calculated as the local optical 
depth towards the identified clump, with kir ~ 5 (Z/Zq) g~' cm^ 
approximately constant over the relevant physical range of dust 
temperatures. The imparted acceleration is directed away from the 
clump, along Rd. There is no "free-streaming" or "tumoff" of cool- 
ing or other physics for particles affected by the feedback. For a de- 
tailed discussion of the normalization, resolution effects, effects of 
photon "leakage," and choice of whether the momentum is applied 
as a continuous acceleration or discrete "kicks," we refer to Paper 
I, but note that essentially all of our conclusions are robust to these 
and other small variations in the model implementation. 

There is also direct momentum injection from supemovae, 
stellar winds, and protostellar jets. We include this in the same 
manner, but with the appropriate momentum injection rates using 
P ~ Pmd + PsNe + Pw , where /'sNc(?, Z) and P„{t, Z) are tabulated 
as a function of time and metallicity for each star particle from 
STARBURST99. We emphasize that this is the direct momentum 
from ejecta only - where SNe are important, they typically impart 
much greater momentum from pressure-driven expansion of hot re- 
gions, which we wish to model directly. The direct ejecta momen- 
tum of SNe and stellar winds are usually small relative to the radia- 
tion pressure (they have a magnitude ~ L/c, but do not receive the 
"boost" of ~ TiR, which is typically at least a few in GMCs). The 
contribution from protostellar jets can in principle also be added 
explicitly, but when integrated over the IMF, their contribution is 
~ 10 — 20% that from stellar winds alone; protostellar jets are prob- 
ably important in regions of low-mass star formation, but they have 
little global effect relative to the contribution of massive stars. 

2.1.2 Supernova Shock-Heating 

In addition to providing momentum in ejecta, the gas shocked by 
supernovae can be heated to high temperatures, generating bubbles 
and filaments of hot gas; the resulting over-pressurized regions can 
then hydrodynamically accelerate nearby gas. The model in § |2.1.1| 
already accounts for the direct momentum of SNe ejecta; but not 
necessarily the full post-shock heating/energy budget (or the mo- 
mentum generated by the pressure-driven expansion of bubbles). 
We therefore add an explicit model for this thermal heating. 

We know the age, mass, and metallicity for each star particle 
in the simulation, so for a given timestep f — )■ / + Af, we can cal- 
culate the total mechanical energy produced by Type-II SNe explo- 
sions. We do so using a pre-tabulated £ (r) (which we compile from 
the STARBURST99 package, for a Kroupa, 2002 IMF; Leitherer] 
|et al.|1999l l. This gives an expected input AZ?. However, on very 
short time/mass scales SNe are discrete events with a roughly fixed 
individual energy, so we actually determine the expected number 
of SNe per star particle in a given timestep assuming a canonical 
Eq = lO^'erg per SNe, and then randomly determine the number 
of SNe in the timestep drawing from a Poisson distribution with 
the calculated expectation value; this ensures the total SNe energy 
comes out correctly. In fact, the time steps in our simulations are 
typically short enough (~ 100 — 1000 yr) to resolve individual SNe 
events. We also include SNe Type-I, for which we adopt an empir- 



ical fit to a simple "prompt-l-delayed" model fitted to the observed 
SNe rates in |Mannucci et al^p006[ l: the prompt component is a 
narrow Gaussian at a few 10' yr, and the delayed component has 
a constant rate. However, since the systems we consider here are 
all star-forming, Type-1 SNe make little difference and the energy 
input is dominated by the Type-II population. 

We wish to explicitly resolve the energy-conserving (Sedov) 
phase of the resulting SNe expansion. We therefore directly deposit 
the SNe energy in the local SPH smoothing kernel (nearest ~ 32 
neighbors), weighted by the kernel - essentially the smallest mean- 
ingful length in the simulation. In this case, there is no implicit or 
sub-grid assumption about even the earliest SNe expansionn 

Unlike most models of supernova feedback, we do not artifi- 
cially turn off gas cooling after injecting thermal energy into the 
ambient gas. Because other feedback processes (radiation pressure, 
photoionization heating, etc.) typically expel the dense gas in star- 
forming regions back out into the more diffuse ISM prior to when 
SNe explode, the effects of SNe on the ISM are significantly larger 
than in models that only include thermal SNe feedback. 

2.1.3 Gas Return and Stellar Wind Shock-Heating 

Stellar evolution returns gas to the ISM from SNe and stellar winds. 
This is implemented following the methodology of |Torrey et al.| 
(j20T2j, with a few minor modifications. Every star has a known 
age, mass, and metallicity; for a given timestep A? we determine 
the total mass-loss rate of the particle (treated as a single stellar 
population) M{t) (tabulated from STARBURST99 for both SNe 
and stellar wind mass losses assuming a Kroupa IMF), which gives 
an expected recycled mass AM = M{t) A?. This mass is then dis- 
tributed among local gas particles in the SPH smoothing kernel, 
weighted by the kernel. Gas recycled in SNe ejecta is returned dis- 
cretely only in SNe (to the same gas heated by those events); gas 
recycled by stellar winds is returned continuously (in every time 
step)|^ Integrated over a Hubble time, the total mass fraction re- 
turned according to the tabulated models is ~ 0.3. 

The injected/recycled gas has the same gross properties (mo- 
mentum, location, etc) as the star from which it spawns. However, 
the initial thermal state must be specified. For the SNe ejecta the 
gas carries the ejecta energy discussed above. We assume that the 
stellar wind gas shocks, so that it has an initial specific energy 



Some other simulations jKatzetal.|l996l [ Stinson et al. |2006| inject the 

energy instead into a larger "blastwave radius" around star particles, with 

p _ <;<: /r0.32_,-0.16 p-0.20 ,,^ 
K£— SDfij, pc (i) 

jChevalier & Gardner|[T974l |McKee & Ostriker|[T9T7) , where £5, = 
Ai?/10" erg, no is the density in cgs, and P04 = 10~^P/[ergcm~^]/:^', 
intended to reflect the radius of the energy-conserving (Sedov) phase of the 
blastwave (often much larger than our softening length). Because we wish 
to resolve this expansion as opposed to assuming it in sub-grid fashion, we 
do not use this approximation. However, we have re-run our simulations 
using this approximation instead of coupling the energy inside of /is,^! . For 
some individual blastwaves it can make a significant difference to their ex- 
pansion history, but on average the behavior is similar where the densities 
are low and bubbles can expand. The final statistical properties of our sim- 
ulations do not depend on the injection radius within a modest range, either 
with this formula or adopting fixed physical radii ~ 1 — lOpc and/or fixed 
neighbor number = 16 — 320. 

' Models in which the gas is stochastically returned by allowing an individ- 
ual star particle to be turned back into a gas particle ("recycled") with prob- 
ability p = l\M/mi give nearly identical results, but are subject to greater 
shot noise, especially on sub-GMC scales. 
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(and corresponding entropy/temperature) set by the wind proper- 
ties as a function of time. Tliis is again tabulated from STAR- 
BURST99 as a function of age and metallicity. Integrated over their 
history, the winds have an energy-flux-weighted average tempera- 
ture r„ « 3 X IO'A", reflecting the fact that most of the energy is 
in early "fast" winds with speeds > 500km s^'; at late times, how- 
ever, the energies and temperatures are much lower because AGB 
winds with velocities ~ lOkms"' dominate. As before, we do not 
modify the cooling or hydrodynamics of these gas particles. 

2.1.4 Photoionization-Heating of HII Regions 

On small scales around young stars, direct photoionization-heating 
by ionizing photons maintains gas in a warm ~ 10* K ionized state, 
which can in principle be an important source of pressure in low- 
density systems. A rigorous solution for the effects of this heating 
would involve on-the-fly radiation transport. This is probably im- 
portant for the study of the formation of the first stars and for de- 
tailed studies of GMC disruption in lower mass molecular clouds, 
but is more detailed than is possible given the goals of the feedback 
model we are developing. Moreover, as we shall show, photoioniza- 
tion feedback is important, but not dominant, for the more massive, 
evolved systems we model. We approximate the key effects of pho- 
toionization on ambient gas with the following simple model. 

Each star particle produces ionizing photons at a rate A', , which 
we again tabulate as a function of the age and mass of the star par- 
ticle from the STARBURST99 models. We then allow these ion- 
izing photons to ionize a region around the star as follows. First 
consider the gas particle nearest the star. If it is already ionized 
(either because it is at sufficiently high temperature > 10"* or be- 
cause it is already a part of another HII region), we move to the 
next nearest gas particle. If it is not ionized, then we calculate the 
rate of ionizing photons needed to maintain the gas as fully ion- 
ized, as AN = N{H) i Prie, v/heve N(H) j =Mj/fimp is the number 
of atoms in the particle, /? ~ 3 x 10^"cm's~' is the recombina- 
tion rate coefficient, and iie is the average election density of the 
particle when fully ionized. If Nj > ANj, the particle is tagged as 
being within an HII region, and the remaining photon production 
rate is tabulated as Ni — >■ A', — ANj. The process is then repeated for 
the next nearest particle. Once we reach a point where A', > and 
Ni < ANj, we determine whether or not to ionize the particle ran- 
domly, with probability p — Ni/ ANj, and consume the remaining 
photon budget, ending the chain. 

This procedure determines the Stromgren sphere around the 
particle, but allows for the resolved inhomogeneous density dis- 
tribution and overlapping ionized regions. In low-density regions, 
the ionized material can extend to very large radii; because we 
also account for a uniform extragalactic ionizing background, we 
also truncate the ionizing region if we reach a radius at which 
the flux density from the ionizing source falls below the ion- 
izing background (for background flux density this is the r 
where L^/47rr~ c = m^). Gas particles flagged as being within an 
HII region are not allowed to cool below an effective temperature 
fffli ~ 10"* K in that timestep (i.e. until they are flagged as outside of 
such a region), and are immediately heated to tmi if their tempera- 
ture lies below this value. 

Although the HII regions estimated for individual star parti- 
cles can be quite small (as the stellar masses can be as low as 
IOOMq), this method allows for nearby star particles to gener- 
ate overlapping HII regions, creating large HII bubbles of sizes 
~ 0.1 — Ikpc, similar to what is observed in many galaxies. Of 
course, if the densities around stars are much higher than what is 



resolved, this model could artificially over-ionize the surrounding 
region; however, given our ability to resolve dense clumps with 
n > lO^cm^^ (inside of which the Stromgren radii are negligible), 
we do not see much resolution dependence until we go to quite 
low resolutions (lower than our "standard" cases). Moreover, a non- 
trivial fraction of photons escape their birth environments, particu- 
larly late in the life of the most massive stars, as they disrupt their 
natal GMCs. 

2.1.5 Long-Range Radiation Acceleration 

Stellar feedback on small scales can disrupt GMCs via the mecha- 
nisms summarized above. On larger scales, galactic-scale outflows 
can be continuously accelerated as well. There are two basic mech- 
anisms that can do this. First, pressure forces from hot gas can drive 
outflows. This should be captured explicitly in our modeling, as 
part of the hydrodynamics, given the heating from photoionization, 
stellar winds, and supemovae. In addition, gas can be accelerated 
outwards by the net radiation pressure produced by all of the stars 
in the galaxy. (Cosmic rays can also contribute to such acceleration 
but are not modeled here.) 

In principle, the physics of radiative acceleration by the cumu- 
lative stellar radiation field is very similar to that associated with 
the local radiation pressure force discussed in § |2.1.1[ but in prac- 
tice, the model in § |2.1.1| reflects the acceleration of gas particles 
by radiation pressure from only the nearest star cluster. These are 
often dense regions in which even the infrared photons are trapped 
so that the local momentum deposition is ~ tirL/c. 

Properly accounting for the full radiative acceleration by all 
of the stars is a formidable calculation, in particular because of the 
possibility of diffusion/scattering in the infrared. Here we make the 
approximation that the dominant acceleration on relatively large 
scales (ie., outside of the dense environments of star clusters and 
their natal GMCs) is due to the optical and UV photons, whose flux 
we can estimate using simple attenuation along the line of sight be- 
tween the gas and star particles. Altogether, the total momentum 
imparted via this absorption/scattering will often be less than the 
momentum imparted locally from the trapped IR photons in the 
star-forming cloud, but the escaping UV/optical photons have the 
ability to accelerate gas at arbitrarily large distances from their ori- 
gin. Moreover, gas above the midplane of the stellar disk "sees" 
the entire surface flux from the disk. Thus the light from many star 
clusters acts together, and can produce a large vertical acceleration, 
in principle launching the gas parcel out of the galaxy entirely. 

For a single stellar source the flux incident on a given SPH 
particle is 

F=-^r (2) 

where r — Tgas — Tsouicc and we have temporarily ignored absorp- 
tion and re-radiation. This radiation flux corresponds to a mo- 
mentum flux per unit area of dP/clAdt = F/c. This is incident 
on a mass m,, area tt/iJ and surface density E, = m,/[7r/i^] (here 
hi — 1.07A'j^^''^/ismi is related to the SPH smoothing kernel, as- 
suming all A'srai particles in the kernel are distributed uniformly 
over the volume and averaging over the kernel mass profile). 
The momentum flux imparted to the gas is just the area times 
the fraction of the flux absorbed or scattered per unit area, i.e. 
P, — {dP/dAdt)Tvhj {i — exp(— kE;)). The acceleration is then 
mi a — Pi,ov 

a=^-V'-— (l-exp(-'^S,)) (3) 
4 IT r' c m, 
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This generalizes to a range of wavelengths using: 

L(l-exp(-KE,)) ^ y^Li,(l-exp(-K^E,))di/ (4) 

If the column through any individual particle is optically thin (typi- 
cal at high resolution), this simply reduces to the standard optically 
thin force term dV/dt = nF /c, and if the column is optically thick, 
it reduces to the assumption that all radiation incident on the sur- 
face 5A = nhj is absorbed. 

In general, seen by a given gas particle is not simply the 
un-absorbed spectrum from each star. It is, however, prohibitive 
to solve the full radiative transfer equations on the fly. The key 
approximation we make is to assume that most of the absorption 
and re-emission happens in the vicinity of each source. The optical 
depth towards a star scales as r oc J pdl ~ E; for realistic density 
profiles the optical depth is often dominated by the contribution 
from small radii near the star, justifying our local approximation. 

In this limit, we can replace in equation [3] with an effec- 
tive L^, which represents the emission of a given source after local 
absorption and reprocessing has been accounted for. The revised 
Li, is modified by local radiative transfer effects on the SED and 
by the fact that reprocessing alters the effective size and geome- 
try of the emitting region. To represent this geometric effect, we 
spread the emission out over a region of length scale with the 
softening kernel (we choose to be the maximum of the local gas 
smoothing or gravitational softening, so it represents both the lo- 
cal gas distribution and implicit stellar mass distribution). We have 
checked that this approximation, in terms of the locations of the 
flux origin, gives good agreement with the results of a full radiative 
transfer calculation including scattering and re-emission performed 
with SUNRISE jJonsson et al.|2009| >. 

The reprocessing and attenuation also effects the emergent 
SED of a given source j via: 

U=Liexp(-A.,E^„,_) (5) 

where LI is the intrinsic (un-obscured) spectrum from the stellar 
population (tabulated as a function of age from the STARBURST99 
models), and Ecoiumn is the total column density (integrated to the 
absorbing particle) from the source. We show in § |Al| that we can 
approximate the total column reasonably well by 

C^LunJ = P.I h.ii ^ Pi [a hj + ^] (6) 

where is the local gas density evaluated at the stellar source, hj 
is the local smoothing length, and the characteristic scale-height of 
the density is given by Vpj|. Note that this second term makes 
the result for Ecoiumn independent of resolution. If the gas density 
around each source fell off with a smooth profile equation|6]would 
be exact, with fi ~ 1 a geometric factor that depends on the density 
profile shape (we adopt ct = 0.95 which corresponds to a linear 
density profile smoothed over the SPH smoothing kernel). 

Because we are only interested in the integrated absorbed mo- 
mentum, it is not necessary to treat as a finely-resolved spec- 
trum. Instead, we integrate the SED over three broad frequency in- 
tervals in the UV (A < 3600 A), optical/near-IR (3600 A< A < 3 ^i), 
and mid/far-IR (A > 3 /i). The "initial" luminosity of a star par- 
ticle Lv is given by the corresponding Luv, i-opt, ^-ir tabulated 
from the STARBURST99 models as a function of age and metal- 
licity. The opacities Hi, are then replaced by the corresponding 
flux-mean opacities kuv, Kopt, kir, which we take to be constant 
= (1800, 180, 5)Z/Zq, respectively|^In general, the opacity will 

In detail, the effective k for each band will depend on the exact SED 



also be a function of the thermodynamic state of the gas, e.g., be- 
cause dust can be destroyed by sputtering in hot gas. These effects 
are not taken into account in our model. 

Given the above model, it is relatively straightforward to im- 
plement the resulting radiation force on all of the gas in the simu- 
lation. Because the acceleration follows an inverse-square law with 
sources of luminosity L instead of mass M, the computation of the 
acceleration is implemented with exactly the same tree structure as 
the existing gravity tree. Of course the force acts only on gas par- 
ticles, with the final normalization coefficient including the terms 
dependent on particle properties (eq.[3](. 

In Appendix|A]we present several tests and calibration studies 
of this model. For example, we show that the approximate treat- 
ment just described recovers the results of full (post-processing) 
radiative transfer calculations of the SED reasonably well. We also 
show that if we were to simply model the radiation field as that em- 
pirically observed in typical galaxies of each type we consider, we 
obtain very similar results for the momentum absorbed. 

We stress that our model of the acceleration by the diffuse ra- 
diation field does not "double count" when combined with the ra- 
diation pressure included locally in § |2.I.I| the radiation pressure 
term there was (I + tir)L/c. The factor "I" comes from the initial 
extinction applied here in equations|5]and[3] exp (— ^column)' i-^-' 
the local absorption of the initial UV/optical photons. In equation 
l[3j, this factor does not add momentum to the gas, rather it atten- 
uates the "initial" spectrum — effectively reducing the momentum 
carried by the UV photons emitted by young stars. It is true that we 
slightly over-count when we approximate the absorbed fraction of 
un-reprocessed radiation by unity in § |2.1.1[ particularly if the ac- 
tual optical depth < 1. However, since the dust opacity 
to the initial radiation (which is mostly in the optical/UV) is large, 
this is rarely the case for the models we consider. In fact, we have 
implemented a more accurate model that allows for UV optically 
thin sight lines in the 'local' (to star clusters) radiation force calcu- 
lation, and find that it produces no measurable differences. 

2.2 Molecular Chemistry 

Another potentially important element in modeling star formation 
is the determination of where gas can or cannot become molecular. 
Stars form only from molecular gas ~ albeit not all molecular gas, 
only that at the very highest densities. We have therefore also im- 
plemented several models which allow for explicit estimates of the 
molecular gas fraction /«, at each point in the simulation. 

In § [5] we compare three different models for the molecular 
fraction and molecular chemistry: 

Molecular Model A: assume /h, = 1 always in regions above 
the SF density threshold. 

Molecular Model B: the model in |Krumholz & Gnedin|pOI l| l, 
which gives a simple fitting function for fn^ as a function of the 
local column density (~ ph^mi) and metallicity. 

Molecular Model C: the model in [Robertson & Kravtsovl 
(j2008}, where both the molecular fraction ///, and the cooling func- 
tion A(r) at temperatures < lO^K are varied as a function of the 



shape, but we choose these frequency intervals in part because the result- 
ing ft are insensitive (changing by < 10%) to the spectral shape for a wide 
range of stellar population ages 1 — lOOOMyr and degrees of dust redden- 
ing. The specific values used here are the average flux mean opacity of a 
solar-composition dust distribution in full radiative transfer calculations of 
the galaxy types modeled here (see AppendixlAl. 
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local density, metallicity, and radiation field to match the results of 
CLOUDY calculations. 

In each case we allow only molecular gas to form stars by 
multiplying the calculated SFR in each particle by . 

We show in §[3] that the choice of models (A)-(C) has no ef- 
fect on our results. This is because we only allow star formation 
in very dense regions (above a threshold n > 100— lOOOcm"^), 
and the galaxies we model here are already substantially metal- 
enriched (even the SMC model has Z ~ 0.1 Zq). Thus the more 
complex models ultimately return the result that the gas in poten- 
tial star-forming regions is overwhelmingly molecular - typical op- 
tical depths to photo-dissociating radiation are ~ 10**, so they are 
strongly self-shielding. 

At sufficiently low metallicities, however, Z < 10^^^ — 
10~^Zq, molecule formation will be inefficient even at high densi- 
ties and cooling times will become comparable to dynamical times. 
The formation of the first stars, and the lowest-mass dwarf galaxies 
(galaxy stellar masses < 10* Mq), are in this regime, where the ex- 
plicit treatment of the molecular fraction above may be important 
( |Kuhlenetal.|2011> . 

3 GALAXY MORPHOLOGIES 
3.1 With Feedback 

Figures [T|5] show the morphologies of each galaxy model, in our 
standard case where all feedback mechanisms are present. 

We show both the gas and stars, face-on and edge-on, at times 
after the systems have reached a quasi-steady state equilibrium. 
The gas maps show the projected gas density (intensity) and tem- 
perature (color, with blue representing cold molecular gas at T < 
lOOOK, pink representing the warm ionized gas at ~ 10* — 10^ K, 
and yellow representing the hot. X-ray emitting gas at > 10* K. 
The stellar maps show a mock three-color observed image, specif- 
ically a u/g/r composite. The stellar luminosity in each band is 
calculated from each star particle according to the STARBURST99 
model given its age, mass, and metallicity (and smoothed over 
the appropriate kernel). We then attenuate the stars following the 
method of Hopkins et al.||2005^ : we calculate the total dust column 
(from the simulated gas) along the line-of-sight to each star particle 
for the chosen viewing angle (assuming a constant dust-to-metals 
ratio, i.e. dust-to-gas equal to the MW value times Z/Zq), and ap- 
ply a MW-like extinction and reddening curve (as tabulated in |Pei| 

The dramatic differences between the different galaxy mod- 
els are obvious. The HiZ models result in disks with massive 
(> 10^ Mq) ~kpc-scale clumps, which form massive star cluster 
complexes. The gas densities are sufficiently high that extinction 
can lead to a clumpy optical morphology as well (by entirely ex- 
tincting out regions in the galaxy). The dumpiness is more promi- 
nent in the gas, but remains visible in the optical images. Seen 
edge-on, the gas can resemble a "clump-chain" morphology. Vi- 
olent outflows arise from throughout the disk, driven by a massive 
starburst with SFR > lOOMQyr^'. 

The Sbc models also produce clumpy disk structure, but the 
individual clumps and star forming regions are relatively much 
smaller, and the global structure more clearly traces spiral arms 
and global gravitational instabilities. As a low-mass starburst (the 
galaxy has a SFR of ~ 1 — 5Mq yr^', despite being an order of 
magnitude less massive than the MW), the irregular gas morphol- 
ogy and thick disk more closely resemble dwarf starbursts in the 
local Universe such as NGC 1569 or 1313. A clumpy, multi-phase 



super-wind is clearly evident, arising from the high specific star 
formation rate. 

The MW models have a familiar spiral and barred spiral mor- 
phologies. But with our high resolution and realistic modeling of 
the ISM, individual star clusters and structure within the spiral 
arms - dust lanes, feathering, and SNe-powered bubbles are evi- 
dent. Edge-on, the dust lanes and moderate/weak outflows of hot 
gas are also evident, but without the violence of the HiZ case. 

The SMC models resemble typical irregular galaxies. The 
ISM is extremely patchy, with distributed molecular clouds and 
large, overlapping SNe bubbles. Edge-on, shells of cold material 
accelerated by winds out of the disk appear as blue "loops" at large 
scale heights, while direct venting of the hot gas is similarly appar- 
ent. The disk is "puffy" and thick, and the stars are distributed in 
the characteristic irregular fashion that defines the morphological 
class. The outer disk is clearly not star-forming; being low-density, 
this is maintained as an extended ionized disk by HII heating (from 
both the inner stellar disk and the UV background). 

3.2 Effects of Each Feedback Mechanism in Turn 

Having seen the typical galaxy morphologies in the case where 
feedback is present, we now consider the galaxy morphologies (and 
star formation histories, discussed later) with different feedback 
mechanisms de-activated in turn. 

Figure |6] shows this for the HiZ and SMC models. We show 
the morphologies (at fixed time and spatial scale) and star formation 
histories of a few models which are otherwise identical but include 
or exclude different feedback mechanisms. Qualitatively, in the HiZ 
model, the clumpy morphology is present in all cases - since the 
galaxy is violently unstable, it is primarily gravity (not feedback) 
that dictates the gas morphology. But there are major differences 
between different feedback models. Most strikingly, without local 
momentum deposition from radiation pressure, the characteristic 
sizes of the individual clump cores are much smaller (this is diffi- 
cult to see by eye - the central regions of each clump that appear 
near- white in the other images are single white pixels in this case!) 
This occurs because the long-range photon momentum, supemovae 
momentum, and stellar wind momentum are not sufficiently strong 
to resist the collapse of clumps beyond a critical threshold. 

The effectiveness of the radiative feedback in the HiZ model 
results from the fact that, as shown in Paper I, the IR optical depths 
reach tir ~ 50. The photon momentum L/c (neglecting the boost 
of Tir), or the direct SNe momentum (also ~ L/c) are lower than 
that produced by the IR radiation by an order of magnitude. Hot gas 
(from the SNe, stellar winds, and HII regions) is clearly insufficient 
to prevent collapse - not surprising since the average density in 
the clumps is > 100cm~^^, so the cooling time is extremely short 
(~ lOOyr). 

Turning off sources of hot gas, the morphologies in the HiZ 
galaxy is similar to the standard, all feedback case - the differences 
being mostly in the diffuse inter-clump gas and diffuse (volume- 
filling) phase of the wind. This gas is much cooler (note the absence 
of yellow gas), and without SNe and stellar winds the clumps and 
spiral arms in the gas appear much smoother (the irregular edge 
morphology and occasional bubbles are largely due to SNe). But 
the absolute wind mass is fairly similar ~ most of it is in the warm- 
phase gas clumps. Heating from HII regions is negligible - unsur- 
prising since ~ lOkms^' is significantly less than the turbulent 
velocities generated by the other sources of feedback. 

We do not show the Sbc and MW cases explicitly (their salient 
properties are shown below), but they are qualitatively similar to the 
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Figure 6. Face-on and edge-on gas morphologies of the HiZ (top) & SMC (hottom) galaxy models with different feedback mechanisms enabled and disabled; 
the stretch and color scale are the same as in Figure[T] Left: Standard model (all feedback on). Middle: No heating: energy injection from SNe, shocked stellar 
winds, and HII photoionization-heating is disabled (gas recycling from these mechanisms remains, but is "cold"). In the HiZ case, the "hot" volume-filling 
ISM is greatly diminished but the global morphology remains similar In the SMC case, the morphology is dramatically different with the heating mechanisms 
turned off - there is no wind and the morphology lacks the characteristic irregular/patchy structure. Right: No radiation pressure (local or long-range). In the 
HiZ case, the hot medium is present, but individual GMC complexes (~ 0.1 — 1 kpc-scale white blobs at left) collapse without limit (to single pixels here) and 
much of the cold wind mass is absent. In the SMC case, GMCs and dense ionized regions also collapse further (leading to the apparently larger filling factor 
of hot gas), but the wind and disk thickness are not strongly affected. 
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HiZ result. In the MW-analogue without radiation pressure, GMCs 
and star-forming "knots" again collapse to significantly smaller 
sizes, although the runaway is not quite as severe (early "fast" stel- 
lar winds and HII heating are able to somewhat suppress runaway 
collapse). Removing the sources of heating does not much change 
the global gas morphology, but eliminates the volume-filling "hot" 
gas and leaves the ISM as a whole much more cold, with a larger 
fraction of material in dense clouds. The She case is lower-mass 
than the MW but more strongly self-gravitating and gas rich, so lies 
as expected somewhere between the MW and HiZ cases in these 
properties. 

However, when we repeat the experiment for the SMC-like 
model (in the bottom half of Figure |6]l, the results are strikingly 
different. The key difference from the other three types of galaxy is 
that the average gas density is much lower, n<0.l cm~^, approach- 
ing the regime where the cooling time is non-negligible even near 
young stars. Moreover, with a lower column density and metallic- 
ity, the IR optical depths can be much lower so photon momentum 
is boosted little beyond ^ L/c. Turning off the sources of photon 
momentum (both local and long-range) has a relatively weak effect 
on the morphology (some GMCs do collapse further, but many are 
disrupted by SNe and stellar winds). 

Turning off all of the "heating terms" from SNe, stellar winds, 
and HII photoionization however, has a large effect - the volume 
filling hot gas "bubbles" disappear and the gas clearly collapses 
along more global gravitational structures into a filamentary mor- 
phology reminiscent of fragmenting very large-ra spiral arms, quite 
different from what is observed in real dwarf galaxies. Considering 
the three heating terms in turn, any one of them is able to prevent 
this complete over-cooling, leading to a morphology more reminis- 
cent of the "standard" case. Quantitatively, however, SNe appear 
to be the most important heating mechanism on large scales: with- 
out SNe, there are still some hot regions created by stellar winds, 
but the integrated energy input in the winds is a factor ~ 8 lower. 
Turning off stellar wind heating has similar but less noticeable ef- 
fects. HII regions and the warm gas pressure they provide are most 
important inside of individual GMCs; they can produce some sup- 
pression of structure but no "hot" gas, and when SNe and stellar 
winds are present, the addition of photoionization-heating has little 
effect on the global gas morphology in the star-forming disk. 



4 STAR FORMATION HISTORIES 

Figure |7] shows the corresponding star formation history for each 
model discussed in § |3.2| with different feedback effects included 
or excluded. For reference, we show both no-feedback models (in 
which case nothing stops catastrophic collapse) and the standard 
all-feedback models for each galaxy type. 

In each galaxy type, with no feedback, the disk experiences 
a runaway collapse leading to the consumption of a large fraction 
of the gas in a single dynamical time (M, ~ Mgas/fdyn). This gives 
much higher peak SFRs (^ 2500, 70, 40, SMqyi^' for the HiZ, 
Sbc, MW, and SMC models, respectively) than observed in the real 
systems the simulations are designed to model, and violates the 
observed Kennicutt relation by factors of ~ 10—100 (see Paper I). 
The SFR declines only owing to gas exhaustion. 

With the standard feedback model enabled, each galaxy type 
reaches a much lower equilibrium SFR where feedback balances 
gravity. These models remain in quasi steady-state for many dy- 
namical times. In Paper I we show these SFRs agree well with 
the observed Kennicutt relation, as is clear from the fact that they 



are in close agreement with the star formation rates observed in 
galaxies of each type, i.e., M, ~ 100 - 300, 1 -5,1-2, 0.01 - 
0.03 Mq yr~ ' , respectively. 

How does this self-regulation depend on the inclusion or ex- 
clusion of different feedback mechanisms? In most cases the differ- 
ent SFHs parallel the differences in morphology discussed in § |3.2| 
For the HiZ model, including or excluding heating from SNe, stel- 
lar winds, and HII regions has a negligible effect on the SFH (re- 
moving all three together is similarly negligible). Excluding either 
of the local momentum or long-range radiation pressure forces in- 
creases the peak SFR by a factor of ~ 2 — 3. Removing both leads 
to a runaway burst similar to the no-feedback case, with a peak SFR 
~ lOOOMQyr"' and then a falloff from gas exhaustion. The sys- 
tem is then an order-of-magnitude above the SFR predicted by the 
observed Kennicutt-Schmidt relation. Clearly the critical physics 
are the cold radiation sources of momentum. 

In the MW-like case, the SFHs are again most strongly af- 
fected by the local momentum deposition and long-range radiation 
pressure, but the effects of SNe heating are non-negligible. Remov- 
ing both local momentum deposition and long-range radiation pres- 
sure forces leads to a near-runaway consumption of gas with a peak 
SFR > IOMq yr~' in the earlier stages of the simulation evolution. 
This is again in excess of the observed Kennicutt-Schmidt relation. 

Removing all nuclear powered sources of heating (SNe, 
shocked stellar winds, and HII photoionization) is somewhat more 
complex: the SFR for the first few dynamical times remains low, but 
then rapidly rises above the "standard" case by a factor ~ 3 — 10. 
The effect is quite non-linear, and in particular non-local. It arises 
from the fact that stellar mass loss is present in the model, but is 
now "cold" (the returned mass is injected at < 10''K and cools 
rapidly). Because the MW is quite gas-poor, this cold recycled gas 
can easily double the available (cold) gas mass supply at small radii 
(the central couple kpc) over the first ~ 0.5 Gyr; what we see is ac- 
tually that the gas lost by stars in the bulge (where the MW is ob- 
served and by construction initialized to have very little gas) builds 
up until a ~kpc-scale gas bar appears and drives a nuclear starburst 
at about ~ 0.5 Gyr. Unlike the "no-feedback" case, the radiation 
pressure terms are sufficient to ensure that the model remains on the 
observed Kennicutt-Schmidt relation at all times: the SFR is high 
because that excess gas becomes compact and builds up a large 
mass, so the system moves along, rather than off, the relation. This 
explains why, in Paper I, with no recycling of stellar gas, we saw 
that the "radiation pressure alone" case was sufficient to maintain 
the MW at a steady-state SFR ~ 1 Mq yr~ ' for several Gyr. 

Among the sources of gas heating, SNe are the most impor- 
tant. If we just turn off shock-heating from stellar winds ("No Stel- 
lar Wind"; recall this retains stellar mass loss but it is cold), the 
equilibrium SFR is higher but only by a factor ~ 1.2 — 2. This is 
not surprising, since the integrated SNe energy deposition is ~ 8 
times larger than that in stellar winds. Turning off HII heating also 
has only a small impact on the equilibrium SFR (when other mech- 
anisms are still present). The HII heating has some effect on the 
evolution of individual star-forming regions (aiding in their disso- 
ciation at early times, if they are relatively low-density), but does 
not strongly affect the global balance of feedback versus collapse. 

As in § |3.2| the Sbc-like case generally falls between the MW 
and HiZ cases. Heating from SNe feedback has a significant effect, 
but weaker than radiation pressure by a factor of ~ 3 or so. 

In the SMC-like case, we see the importance of hot gas re- 
flected in the SFHs. Removing HII heating has a negligible ef- 
fect on the SFH, and removing hot stellar winds increases the 
peak SFR by a small factor. Removing both the local momentum 
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Figure 7. Stai' formation histories (total SFR versus time) for each of the galaxy models discussed in § |2] with no feedback (dotted line), all feedback 
mechanisms enabled (the "standard" model, solid line), and various mechanisms turned on and off in turn. In every galaxy type, with no feedback the disk 
collapses into GMCs that turn into stars in runaway fashion, leading to enormous SFRs exceeding the Kennicutt-Schmidt relation by ^ 10— 100. These 
enormous SFRs decline over time, but only because most of the gas is consumed into stars. By contrast, with all feedback enabled, the disks reach an 
equilibrium SFR and maintain a quasi steady, feedback-regulated equilibrium (the star formation still declines with time in the HiZ and She models because of 
gas exhaustion). Top Left: HiZ disk. Radiation pressure dominates the regulation of the SFR; turning off all heating mechanisms (dashed line) makes a small 
~ 20% impact on the SFR. Top Right: Sbc starburst disk. Radiation pressure is still very important, but hot gas pressure makes a significant contribution to 
self-regulation. Removing all the explicit heating mechanisms (dashed line) leads to factor ~ 3 higher equilibrium SFR; removing radiation pressure (dot- 
dashed line) leads to factor ~ 10 higher SFR. Bottom Left: MW-like disk. The interplay between the different feedback mechanisins is more subtle in this case 
- see the main text (Q for a detailed discussion. Bottom Right: SMC-like isolated dwarf. Here, SNe heating is the inost important mechanism regulating the 
galaxy-scale star formation. 



and long-range radiation pressure forces only increases the median 
SFR by a significant but not dramatic factor ~ 2 at early times. In 
contrast, removing SNe increases the typical SFR by a factor of 
~ 5 — 8. It is worth noting that even after this increase, the SFR is 
still much less than it would be with no feedback (~ 0. 1 Mq yr^ ' 
versus ~ 3MQyr~'), so there are important contributions to self- 
regulation from all sources of feedback, but SNe are the most im- 
portant among them. 



5 THE (NEGLIGIBLE) ROLE OF MOLECULAR 
CHEMISTRY 

It is well-established that stars form only from molecular gas 
- albeit not all molecular gas, but only a small fraction at the 
highest densities. Hence a number of groups have attempted to 
model galactic star formation using various schemes to estimate 
the molecular fraction in a large volume of gas in a simulation (see 
IRobertson & Kravtsov||2008l |Ostriker et aL||20T0| |Kuhlen et al.| 



|20I \\ . Given our inclusion of a wide range of feedback physics, 
does chemistry matter? 

We test this by re-running our standard models with all feed- 
back present, but with different explicit models for the molecular 
fraction in gas. First (Molecular Model A), we make the crude as- 
sumption that the molecular fraction is always unity in all regions 
above our density threshold for star formation. 

Second (Molecular Model B), we adopt the prescription in 
[Krumholz & Gnedin|pOTT| eqn. 1-4 therein), which approximates 
the results of full radiative transfer calculations with a simple an- 
alytic fitting function that gives (the fraction of gas mass in a 
given particle which is molecular) as a function of the local column 
density and metallicity. We approximate the local column as the in- 
tegral across the SPH smoothing kernel (oc ph, analogous to what 
is done with the AMR experiments therein)]^ We then allow only 

^ This is a conservative assumption because, if the integral column over 
large spatial scales is much larger, then it will further self-shield the gas, 
which only strengthens our conclusions below. 
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Figure 8. Star formation history for SMC models with different treatments 
of the molecular chemistry of the gas; all of the simulations have all of our 
feedback mechanisms enabled. The SMC model is the most sensitive of 
our galaxy models to the molecular chemistry because of its lower surface 
density and metallicity. (A): We use our standard cooling curve (see ^2} 
and assume that all gas above the star- forming threshold (n 2> lOOcm" ) 
is molecular (the H2 molecular fraction /«, = 1) and so can turn into stars 
at an assumed efficiency per dynamical time of 1.5%. (B): We use the fit- 
ting functions in |Krumholz & Gnedin|j201 1| for ///j as a function of local 
column density and metallicity for each gas particle and only allow star for- 
mation in molecular gas. (C): We use the model in |Robertson & Kravtsov| 
j2008| , in which the molecular fraction and the cooling function A(T) 
at low temperatures < lO'* K are a function of the local gas density, metal- 
licity, and radiation field (based on CLOUDY calculations); star formation 
again only occurs in molecular gas. Overall, the treatment of the molecular 
chemistry has a negligible effect on the star formation history. This is be- 
cause star formation only occurs in high density gas, n> 100 — lOOOcm^', 
and at the SMC's metallicity of Z ~ 0. IZq, the molecular fractions in 
the dense, star-forming gas are always near unity (average //^j > 0.998 for 
models (B) and (C)). In test calculations, only at extremely low metallicity, 
Z <C IO~^Z0, does the molecular' chemistry begin to have a measurable 
effect on the star formation history. 



a large volume of the ISM that is being estimated ~ radii at least 
> 100 pc and as large as several kpc, with bulk volume densities of 
~ 0.1 — 10cm~^. Essentially, then, the "molecular" model is really 
an implicit model for the fraction of dense gas within the patch. In 
contrast, our simulations explicitly resolve dense regions, and only 
allow star formation in them; our typical threshold for star forma- 
tion is a local density of n > 100 — lOOOcm"^. At these densities, 
the gas is overwhelmingly molecular, so explicitly accounting for 
a molecular fraction and allowing only that gas to form stars is es- 
sentially identical to our existing prescription. In fact, we know the 
typical optical depths through the star-forming patches of gas (from 
our local momentum-driving model) and in the optical and UV 
wavelengths, they are typically ~ lO"*, so there is no question that 
the gas can self-shield. In fact, when we implement the model of 
IKrumholz & GnedinlpOTT) , we find the average estimated molec- 
ular fraction in the gas we would allow to form stars is > 99.8%. 

Furthermore, in Paper I, we show explicitly that the global 
star formation history is independent of the small-scale star forma- 
tion law, because the star formation rate is set by the requirement 
that feedback from young stars maintain marginal stability. This 
requirement of vertical pressure support implies that allowing only 
molecular gas to form stars on a constant timescale within these 
high density regions will give the same star formation rate. 

We also have tested, in Appendix A of Paper I, that chang- 
ing the cooling curve shape dramatically - for example, changing 
its normalization by a factor of several, or changing its shape to a 
simple constant from 10— 10"* K - has no impact on our results, 
because in each case the cooling time in the gravitationally col- 
lapsing dense regions is still shorter than the dynamical time by a 
factor ~ 10* — 10^ at densities near the threshold for star formation 
(at least down to ~ lOOK, at which point the thermal pressures are 
totally negligible compared to gravity, and so subsequent cooling 
has only weak effects). 



the molecular gas to form stars - multiplying the SFR our standard 
code would estimate by /h^ . 

Third (Molecular Model C), we adopt a more detailed treat- 
ment motivated by that in |Robertson & Kravtsov|j2008^ . We allow 
the cooling function A to be a function not just of temperature but 
also the local gas density, metallicity, and radiation field, where the 
latter is approximated by the sum of the cosmic background (z~0 
value) and the local SFR (taking the instantaneous SFR in the local 
gas and assuming the continuous SF limit) with the SEDs tabulated 
in STARBURST99. Specifically the cooling functions are tabulated 
from Figure 1 of [Robertson & Krav tsov (2008) (on a simple grid 
in density, metallicity, and radiation field) and we interpolate loga- 
rithmically. We then calculate the local molecular fraction in each 
particle as a function of column density, metallicity, and also ioniz- 
ing field strength from the fitting function in Krumholz & Gnedm] 
(j20TTJ. As before, we then only allow SF from the "molecular" gas. 

Figure [8] shows that modeling the molecular gas fractions has 
no effect on our results. If we plot the morphology of each system 
(or any other measurable quantity presented in this paper), we find 
the same. 

Why do we find that detailed chemistry is essentially irrele- 
vant, while other authors have found it is critical? The key is the 
scales that we resolve and on which we allow star formation to pro- 
ceed. In previous models where the molecular fraction is estimated 
and makes a large difference, it is usually the molecular fraction in 



|Glover & Clark| ( |2012[ l study star formation on small scales 
with a variety of different chemical and cooling models, and reach 
the same conclusion we do here. They show that the molecular frac- 
tion under otherwise identical densities and metallicities is causally 
irrelevant to the collapse of dense gas, since cooling via fine- 
structure and dust emission is nearly as efficient as CO - the pres- 
ence of CO is very tightly correlated with, but does not drive the 
collapse of gas to high densities. 

Thus it is only if cooling and the formation of metals, dust, 
and molecules are dramatically suppressed that the explicit chem- 
istry will become important for any of the quantities we model. 
This occurs at metallicities below Z ^ 10^' — 10~" Zq. Specifi- 
cally, this is the metallicity at which the molecular fraction begins 
to fall below unity in gas near our threshold density; this is also the 
metallicity at which the cooling time can become comparable to, 
or longer than, the dynamical time (see e.g. [Robertson & Kravtsov| 
|2008[ l. The formation of the first stars, and the lowest-mass dwarf 
galaxies (galaxy stellar masses < IO^Mq) are within this regime 
and will be strongly influenced by the ability of the gas to form 
molecules even at high densities. But the more massive systems we 
model here (even the SMC) are safely in the regime in which the 
dense gas is almost all molecular and thus so long as star formation 
is restricted to the dense gas, an explicit treatment of the molecular 
chemistry is not important. 
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6 GLOBAL PROPERTIES: DISK VELOCITY 
DISPERSION AND SUPPORT 

We now examine the structural properties of our model disks, as 
a function of time and radius, with different feedback mechanisms 
enabled and disabled. 

Figure|9]plots a number of disk structural properties as a func- 
tion of time. We are particularly interested in the properties of the 
star-forming disk, as opposed to e.g. the extended ionized HI disk 
(which is obviously not strongly affected by feedback), so we cal- 
culate all values as gas-mass weighted averages in each radial an- 
nulus at each time, then plot the SFR-weighted average over ra- 
dial annuli (the results reflect the value near ~ Rc, but do not vary 
rapidly in radius). First, the resolved vertical velocity dispersion 
a- (specifically, half the 15 — 86% interval in V;, to minimize the 
bias from large "tails" of outflowing material at high-V;). Second, 
the ratio on/a-, between the in-disk plane radial velocity disper- 
sion and vertical dispersion (the azimuthal dispersion ~ gr in 
all cases). Third, the (fractional) disk scale height h/R, where the 
scale-length h of the disk fit to a vertical sech^ profile. Fourth, the 
Toomre Q parameter. 



5vrk 
ttGE 



(7) 



where 5vr is the effective radial dispersion of the gas, = 



; (jR (because the thermal sound speed c, <C ctr in the 



dense gas, in all cases), k is the epicyclic frequency (~ ^/iQ, 
for nearly constant-V^ disks here), and E is the gas-l-stellar disk 
surface density. Finally, the disk gas fraction, defined locally as 

/gas = Egas/(Egas + E«). 

The very early time behavior is a result of the initial con- 
ditions: these include large vertical support for the disks, which 
quickly vanishes as the disks cool, then dense GMCs form and be- 
gin to produce stars. After about an orbital time, the systems reach 
a quasi-equilibrium state, which can clearly be seen (at least in 
cases with feedback) in the flat or slow subsequent time depen- 
dence. Once in this regime, details of the initial conditions (e.g. 
whether the initial disks are supported by turbulent or thermal ve- 
locities) have no significant effects. 



clumpy and self-gravitating; feedback and gravity therefore main- 
tain large dispersions in the gas. 

The MW-like runs maintain the canonical observed ~ 8 — 
15kms^' dispersion, which with V^ax ~ 220kms^' gives V/ct ~ 
20 — 30, and gas fractions ~ 10% - because they are much more 
stable (with large dark matter halos), and have low gas fractions 
(and low specific star formation rate), these disks are thinner, rela- 
tive to the starburst systems. 

The SMC-like runs show a similar moderate dispersion ~ 
6 — 8kms^', giving V/a ~ 5 — 15, and maintain very high gas 
fractions ~ 70 — 80%; these properties all agree very w ell with 



typical properties for ~ 10 Mq dwarf galaxies (see e.g. Bell & 
|de Jong|2001[[McGaugh|2005[ l. Here, the self-gravity, star forma- 
tion rate, and corresponding dispersion are low, but the potential is 
very shallow, so even small dispersions maintain relatively "puffy" 
gas disks, as observed. Gas fractions are generally larger at large 
radii as the disks have some bulge and a more extended gas versus 
stellar disk; there are some small flares at large radii but these are 
weak. There is a rise in V ja at large radii as the gas at these radii is 
very low-density and has a low specific star formation rate, hence 
low feedback support and low a. In these outer regions the disk 
is supported in part supported by thermal pressure, maintained by 
external radiation, either from the interior of the disk or from the 
cosmic background ultraviolet radiation. 

We have also measured the full turbulent velocity power spec- 
trum in these simulations, and find its shape - essentially consistent 
with equal power on all scales (with a break below the disk scale- 
height where the system transitions from being effectively two- 
dimensional to three-dimensional) - varies little across simulations 
(with normalization just set by a and the scale-lengths here). This 
is expected: since the bulk gas motions are super-sonic and gravity 
is the dominant force, the turbulence should be largely scale-free. 
This is true even if motions are "pumped" by stellar feedback, be- 
cause the fraction of the time that direct feedback forces are dom- 
inant over gravity and bulk hydrodynamic forces as the force on 
a given gas parcel is small. For this reason, we do not see strong 
imprints of the input feedback injection or velocity scale - quite 
unlike a case where gas is "held up" in some hydrostatic sense by 
feedback forces. 



6.1 Dispersions and Gas Fractions 

6.1.1 Comparison to Observations 

For the cases with all feedback mechanisms enabled (our "stan- 
dard" model), the values generally agree with observed disks at 
similar mass and redshift to the analogs in each model. The HiZ 
disks have SFRs of ~ 70 — 300 Mq yr^ ' , a morphology of massive 
clumps, purely vertical velocity dispersions of ~ 30 — 40 km s^' 
and sightline-averaged dispersions of ~ 50 — lOOkms"', giving 
typical 1//ct ~ 5 or so given their maximum circular velocity 
Vmax ~ 230km s~', and high gas fractions ~ 30 — 60% for most 
of their lifetimes|^all in good agreement with observed properties 
of the most massive star-forming disk s jGenzel et al.|2008[[F6rster| 
ISchreiber et al.|2006[|Erb et al.|2006||Tacconi et al.|2010) . 

The Sbc disks also maintain high gas fractions, low V /o ^ 
5 — 8, with large dispersions and puffy scale heights. In both of 
these cases, the specific star formation rate is high, and the disks are 

* Recall this is essentially /gas at the effective radius of star formation, 
which evolves weakly in time as gas can be depleted at small radii (giving 
the slow decline in SFR seen in FigurelTl. 



6.1.2 Failings of Feedback-Free Models 

There are important differences if we compare runs with feedback 
to the no-feedback cases. Most obviously, without feedback the 
gas fractions plummet extremely quickly, as essentially the entire 
galaxy gas supply collapses into dense regions and is turned into 
stars in a single global dynamical time. 

The vertical dispersions (and as a consequence, scale heights) 
are also significantly lower (at least for most of each run) in the 
cases without feedback. This is despite the fact that the rapid gas 
depletion should make it much easier for non-feedback (and non- 
nuclear powered) mechanisms to maintain large dispersions in the 
remaining gas. In the HiZ and Sbc cases, the systems are strongly 
gravitationally unstable, with large clumps scattering off each other 
and massive inflows along bars and spiral arms. However, almost 
all of this gravitational potential energy goes into maintaining the 
in-plane dispersions. The in-plane dispersions can be maintained 
at large values by global gravitational instabilities (sinking clumps, 
bars, and spiral arms) that avoid rapid shocks and dissipation pre- 
cisely because they are global structures. In contrast, the vertical 
motions drive shock heating on the local dynamical time, so o- 
rapidly drops to a low value. 
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Figure 9. Global galaxy structural properties in gas (not stars) as a function of time, in otherwise identical simulations of each galaxy model with various 
feedback effects enabled (different lines - see Table 2). Each quantity is the weighted average over the star-forming disk. Top: Sightline-averaged one- 
diinensional gas velocity dispersion cr. Second from Top: Velocity (an)isotropy: ratio of in-plane radial dispersion to vertical dispersion cr-. Middle: Vertical 
scale height-to-radius ratio ft/R. Second from Bottom: Toomre Q in the gas (o-rk/ttCB). Bottom: Gas fraction Mgas/ (Mgas +M»). With all feedback enabled, 
values agree well with those observed: /gas evolves slowly, feedback preserves stability (Q ~ 1, with a driven to the appropriate observed value for each 
galaxy) and isotropizes the dispersions (ctr a^, giving h/R a^/Vc). With no feedback, runaway star formation consumes the gas. Dispersion and Q still 
increase, but only because all high-density (rapidly cooling, Q < 1) gas is totally exhausted. In the HiZ and She cases, removing radiation pressure yields 
results similar to the "no feedback" case (SNe, stellar winds, and photoionization-heating cannot prevent runaway dissipation in dense gas). In the SMC case, 
gas heating by SNe is more important than radiation pressure for the disk structure. In the MW case, all mechanisms contribute comparably. 
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Figure 10. Face-on map of the Toomre Q parameter (eqn. |7j for the gas 
within the disk averaged over ~ lOpc annuli. Intensity is scaled logarith- 
mically with Q from Q < 0.1 (black) to 2 > 10 (white). We show the SMC 
model as Fig.|5] Despite the robust global convergence to Q ~ 1 in Fig. [9] 
there are large local variations in Q, as some regions have dissipated their 
support and are collapsing (on their way to forming stars) while others have 
just expeiienced strong local bursts of stellar feedback and are rapidly ex- 
panded and/or heated. 

At late times in both the HiZ and She cases a. begins to rise, 
but the models have a have a gas fraction of just a few percent at 
the time, and that gas is mostly being "dragged" by a few super- 
massive stellar clumps, and undergoing mergers at the galaxy cen- 
ter. 

The SMC case is physically quite different from the HiZ and 
She models: the system requires only small dispersions, but it is 
so gravitationally stable (owing to the large dark matter to disk ra- 
tio) that it does not develop large turbulent motions in the absence 
of feedback. Again, though, it shows a large cr^/cr, without feed- 
back, from the moderate spiral structure in place. The MW model, 
interestingly, is in some ways a "worst-case" system for discrimi- 
nating the effects of feedback (the differences are minimized) - this 
is because it has a low gas fraction, requires a low gas dispersion 
for vertical pressure equilibrium, and is marginally unstable to e.g. 
spiral arm formation. 



6.1.3 Effects of Each Feedback Mechanism In Turn 

Now consider cases with different feedback mechanisms included 
or excluded. For the HiZ and Sbc-like cases, removing the sources 
of "hot" gas (e.g. SNe feedback, stellar wind shock-heating, and 
HIl heating) has relatively little effect on any of the quantities plot- 
ted in Figure |9] Removing the radiation pressure (momentum feed- 
back acting on cold, dense gas) however, leads to much more rapid 
star formation and hence a runaway decrease in /gas. The full dis- 
persion, (j^ + + o-j , is essentially the same for any combina- 
tion, or lack, of feedback. As we argued in the previous paragraph, 
global gravitational instabilities (powered by accretion through the 
disk) can maintain and gr. But clearly the no-feedback and ra- 
diative momentum models are "failed" in the sense that they did 
not prevent runaway star formation and, at a fixed gas fraction fg, 
the vertical velocity dispersion a- is too small. 



For the SMC case, "hot" feedback is much more important, 
since the average densities are much lower and opacities and spe- 
cific star formation rates (important for radiation pressure) much 
lower as well. Removing all the hot gas sources, but keeping ra- 
diation pressure in place, moderate but slightly lower vertical dis- 
persions ~ 5 — 6 km s ^ ' can be maintained, and star formation effi- 
ciencies are still much lower than the no-feedback case (by a factor 
--^10, see Paper I and Q but they are higher than the case with hot 
gas sources by a factor of several. Supemovae are most important 
source of hot gas in the SMC model (as we saw for the SFR above 
and winds in |Hopkins et al.|201 Ib^ . Removing stellar winds or HII 
heating has very little effect on the dispersions or other properties 
if SNe feedback is still included; however, if SNe are absent, they 
are able to maintain comparable dispersions. 

Finally the properties of the MW-like model lie in between the 
HiZ and She models on the one hand, and the SMC model on the 
other, as expected. As noted above, in the MW model all the feed- 
back mechanisms seem to contribute comparably, making it diffi- 
cult to distinguish the effects of different feedback mechanisms. 



6.2 Toomre Q and Disk Stability 

All of the systems modeled here equilibrate at a Toomre 2 ~ 1, and 
maintain this over the disk surface. In the feedback-regulated cases, 
this must happen: Q < \ regions will collapse and begin forming 
stars, which leads to a super-linear feedback (in that the further col- 
lapse proceeds, the more rapid the SF and hence momentum/energy 
injection becomes), until it arrests further collapse (which by defi- 
nition is when it re-equilibrates to ~ 1). If 2 > 1, then there is no 
collapse, hence no dense regions form and no star formation takes 
place; without feedback to hold up the disk, the cooling time is 
much less than the dynamical time, and so it radiates away its sup- 
port in at most a crossing time (the dissipation time for turbulent 
support). 

Fig. [To] illustrates how this global equilibrium is maintained, 
by showing a map of the local values of Q (calculated within lOpc 
annuli) for the SMC model (the other models are qualitatively sim- 
ilar). Although the disk maintains an average 2 ~ 1 there are large 
local variations from 2 ~ 0.1 — 10, corresponding at low-2 to the 
regions which have begun to collapse, and at high-Q to regions 
which have recently formed massive stars and are being disrupted 
and heated by feedback. The characteristic spatial and time scales, 
therefore, for these fluctuations, are closely associated with the 
Jeans scale and associated crossing times. 

However, even the feedback-free systems also equilibrate at 
2 ~ 1. The reason is simple: nothing can prevent runaway collapse 
and star formation in 2 < 1 regions, so all the gas that can reach 
2 < 1 does so, and then turns into stars, leaving only 2 si 1 g^s re- 
maining. Thus 2 ~ 1 can be maintained, but at the expense of run- 
away gas consumption. Quantitatively, we make a crude analytic 
estimate of the dispersions that can be maintained by inflow and in- 
stabilities. The vertical velocity dispersion dissipates on a vertical 
crossing time, H /cr^ = r/Vc, where the second expression follows 
from the assumption that the turbulence provides the vertical sup- 
port against gravity required to maintain the vertical scale height H. 
The luminosity dissipated by turbulence is then L,urb = rja'^Q.Mg, 
where Mg is the gas mass of the disk and 77 is a constant of or- 
der unity. In the absence of feedback, this luminosity must be sup- 
plied by global gravitational instabilities (sinking clumps, bars, spi- 
ral arms), ultimately powered by accretion of gas through the disk. 
(We assume that the accretion energy is converted into isotropic 
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turbulence). The accretion luminosity is Lace = 2Moc.c.V^^.|^We can 
find a lower limit to Mace assuming the disk is in a quasi-steady 
state, so that the mass accretion rate must exceed the star forma- 
tion rate, given by the Kennicutt-Schmidt relation, E, = eEgasfi, 
with e — 0.017. Averaging over the disk, M, = €U~^MgQ where Q, 
is evaluated at the disk effective radius, and ~ 2 depends on the 
surface density profile. Combining these relations, we find 

\2r,uJ 



a- 



(8) 



Using this estimate for a- in the expression for the Toomre param- 
eter. 



2k 



Vc > vh. 



1k \ / 3e yn 

uQJ \2ur]/ 



/gas 



(9) 



so we expect 2 ~ 1 as /gas depletes below < 0.1, which is roughly 
what the simulations show. 

Even in simulations with no star formation, and hence no 
gas consumption ("No Feedback or SF" in Table |2]l, the systems 
quickly reach 2 ~ 1. In these models, the gas collapses and cools 
until it is essentially all in extremely dense clumps, n > 10*cm~^ 
in Fig. 3 of Paper I - corresponding to ~kpc-scale patches col- 
lapsing to ~pc. This makes the gas effectively collisionless! The 
gas "nuggets" (which internally continue to dissipate, but are long- 
lived and rarely collide) scatter and develop Q ~ 1 just as the stars 
do. 



low densities - and although it contains some cold, dense clumps 
(visible in the density distribution or in Figures [T|5) , it contains 
most of the hot gas with Tcff -SJ 10* AT. Much of this material is out 
at very large radii of order the virial radius. 

We next consider the phase structure of the star-forming 
disk in detail, since this is the traditional region in which the 
ISM phases are defined and where stellar feedback has the great- 
est impact. Specifically, we show the mass-weighted density dis- 
tribution (dm/dlogn) and volume- weighted density distribution 
(dV /dlogn; where the volume of each particle is proportional to 
its smoothing length cubed), within the star-forming disk identified 
above (since there is no outer boundary to the simulation, the "to- 
tal" volume-weighted density distribution diverges at low densities 
with the material at the furthest edges of the simulation). We di- 
vide the gas into three traditional phases with a simple temperature 
cut: "cold atomic-l-molecular" gas (T < 2000 A'), "warm ionized" 
gas (2000 < r < 4 X lO^AT) and "hot" gas (T > 4 x K). The 
exact temperature cuts are arbitrary but the qualitative comparison 
does not change if we shift them by moderate amounts. Unsurpris- 
ingly, the high-density gas (> lcm~^) is predominantly cold, the 
intermediate-density gas (0.01 — 1 cm~^) primarily warm, and the 
low-density gas (< 0.01 cm~') hot. The temperature structure also 
explains the features seen in the total density distribution, arising 
from the separate components ~ the total distribution is crudely log- 
normal, but each component individually is more clearly so. 



7 ISM PHASE STRUCTURE 

7.1 Mass and Volume of Phases with Feedback 

Figure pT] shows the density distribution in three of our "standard" 
simulations with all feedback implemented normally. For each, we 
present the distribution in three ways. First, we consider the total 
mass-weighted density distribution (dnigas/ dlogn), for all gas in 
the simulation. This covers a very wide dynamic range in both den- 
sity and spatial scale, and is not directly comparable to the typical 
"ISM density distribution" from galaxy studies. To see this, we di- 
vide the gas into three categories: the star-forming disk (defined as 
gas within the projected radius that encloses 90% of the total star 
formation and within one vertical scale height of the disk inside that 
radius), the extended disk-l-halo gas (gas outside the disk, but not in 
the wind defined below), and the wind/outflow material (gas above 
a scale height, with large radial outflow velocity Vr/|v| > 0.75 and 
Vr > lOOkms"'). The star-forming disk (radii ~ 8, 4, 10, 3kpc in 
the HiZ, She, MW, and SMC case) extends out to a couple opti- 
cal scale-lengths, and (unsurprisingly) includes almost all of the 
dense gas (n > 1 cm~^^). The extended disk and halo cover interme- 
diate densities: this includes the initial extended gas disks, which 
extend to very low densities and tend to be ionized hydrogen disks 
with Tetf ~ 10^ if, and material which is pushed out of the disk by 
feedback but does not contain enough energy/momentum to drive 
an outflow (the "halo" - recall there is no initial extended gaseous 
halo in these simulations). This contains a large fraction of the gas 
mass in all cases, most of which is out at several optical disk scale 
lengths. The wind material extends to the lowest densities - since 
there is no IGM, in principle escaped material can reach arbitrarily 

^ The factor of four (2 = 4x 1/2) comes from the energy transported by the 
torques causing the accretion; this factor is larger (by 4/3) in a constant-Vr 
disk than that in a Keplerian disk. 



By mass, the cold/molecular phase dominates, but with a com- 
parable contribution from the warm/ionized phase, and the ratio of 
the warm-to-cold gas increases as we move to lower-mass, low- 
average density systems (M„anii/Mcoid ~ 0.3, 0.5, 0.7 in the HiZ, 
MW, and SMC-like cases, respectively). The "hot" phase consti- 
tutes at most a few percent of the total mass budget within the star- 
forming disk (~ 1 — 3% in each case). By volume, however, the 
results are reversed. The "hot" phase has a near-unity volume fill- 
ing factor in each case. The "warm" phase has a non-trivial filling 
factor (~ 0.3, 0.5, 0.7 in the HiZ, MW, and SMC-like cases), but it 
can be significantly less than unity. Recall from Figures [T|5] that the 
"warm" phase is still often concentrated in dense structures such 
as the spiral arms (or in the HiZ case, in giant cloud complexes) 
as well as loops and shells blown by hot bubbles (although again, 
most of the warm gas is in the extended ionized disk). The cold 
phase has a very small filling factor of a few percent (~ 3 — 10%). 

The warm and hot phases of the ISM have similar thermal 
pressure P ~ 5 x 10~'^ erg cm""', which for the warm material re- 
flects gas with « ~ 0. 1 — 0.4 and r~l— 5x10"* and for the hot gas 
n ~ 0.003 with T ^ 10* K. Both include a tail of material with pres- 
sures up to ~ 100 times larger, much of which is recently-heated 
material that is venting out of the disk, but some of which is re- 
cently photo-ionized dense gas that is "breaking out" of GMCs. As 
expected, the thermal pressure of the cold gas is lower than the am- 
bient medium by a factor of ~ 10 — 30. However, the turbulent pres- 
sure of the clouds is much larger, P ~ 10^'" — 10"' erg cm"' - this 
is consistent with their being marginally bound (discussed below), 
for which P ~ GE,-,„,d ~ lO-'ergcm"' (Edoud/100MQpc-^)l 
The warm/hot pressure increases as well if we include the turbulent 
ram pressure, but by order-unity factors; they satisfy 2 ~ 1 hence 
P ~ G{T,y, so the pressure ratio is approximately the squared ra- 
tio of GMC to disk- average surface densities. The cold clouds are 
therefore clearly not typically pressure-confined, but held together 
by self-gravity. 
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Figure 11. Density distribution of different ISM pliases, for different galaxy disk models (HiZ, MW, and SMC), at times after they reach feedback-regulated 
steady-state. Variation in time after the first orbital time is modest. Top: Mass-weighted density PDF (dmgas/dlogn), i.e. the mass fraction per logarithmic 
interval in density n. We show the distribution for all gas in the simulation (black), the gas within the (multi-phase) star-forming disk (blue), the gas in the 
extended, ionized non-star forming HI disk and halo (green), and wind/outflow material (orange). These trace the material at low, intermediate, and high 
densities, respectively (as expected). Each density PDF has a very broad density distribution. Middle: Mass- weighted density PDF within the star-forming disk 
- i.e. the traditional multi-phase ISM. We show the density PDF for all of the gas in the disk (black), the "cold/molecular" phase (purple; T < 2000 K), the 
"warm ionized" phase (cyan; 2000 < T < 4 X 10^ K), and the "hot diffuse" phase (red; T > 4 X 10^ K). Most of the mass is in the cold phase (which dominates 
at high densities), but with a non-negligible (~ 20 — 30%) contribution from the warm medium (the hot phase contributing ~ 5 — 10%). There is also a weak 
trend towards the cold phase being more important in more massive, gas-rich disks (HiZ versus SMC). Note that the multi-modal nature of the total density 
PDF is a consequence of the strong phase separation. Bottom: Volume-weighted density PDF (dVgas/dlogn) for the star forming disk. The hot diffuse phase 
dominates, with a moderate (~ 20 — 40%) volume filling fraction for the warm phase, and a small (~ 1 — 5%) volume filling fraction of cold clouds. 



7.2 Dependence on Feedback Mechanisms 

It is useful to identify observational tests which discriminate be- 
tween different feedback mechanisms, and in this section we ex- 
plore some candidates. In Figure[T2] we consider how the distribu- 
tion of cold and hot gas depends on the various feedback mecha- 
nisms implemented. With no feedback, there is runaway collapse 
in the cold gas with a large secondary peak at n — >■ 10^cm~^ (the 
highest densities we can resolve). Turning off all heating mech- 
anisms (SNe, stellar winds, and HII photoionization-heating) to- 
gether leads to a modest increase in the amount of high-density 
n > lO^'cm"^ material in the MW and SMC models (these mecha- 
nisms are also responsible for heating the small tail of "cold" ma- 
terial at n ~ 0.01 — 1 cm~' which appears when feedback is turned 
off). But in all models, turning off radiation pressure yields a much 
more dramatic increase in the amount of very dense material. In 
other words, even where global self-regulation can be set by SNe 



heating, the dense material at n > 10''cm^^^ is regulated by radia- 
tion pressure. This should not be surprising - even in a low-density 
galaxy like the SMC, the optical depths at these densities are large 
in the infrared; and the cooling times for a SNe remnant in such 
high-density material are ~ 10"* times shorter than the dynamical 
time, so pure gas heating is ineffective. 

With no feedback, there is essentially no hot gas, in any model 
- pure gravitational turbulence and shocks alone can maintain some 
gas at ~ 10" K or so, but not > 10* K. Without gas "heating" from 
SNe and shocked stellar winds the mass of hot gas drops by a large 
factor of ~ 10 — 100, even in the HiZ and Sbc cases. The turbu- 
lent motions in the disk driven by radiation pressure have velocities 
comparable to the escape velocity from massive star clusters - tens 
of kms^', so shocks between them produce little "hof material. 
HII photoionization-heating obviously does not produce tempera- 
tures in this range. As we have shown above, most of the hot gas 
comes from SNe (the integrated hot gas energy in stellar winds 
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Figure 12. Top: Mass-weighted density distribution of "cold" gas {T < 1000K)for models with different feedback enabled (as in Fig.jvjand Table 2). Bottom: 
Mass-weighted density distribution for the "hot" gas (T > 10^ K). Without feedback, gas piles up at > 10^ cm~^, and there is essentially no hot gas. Removing 
radiation pressure leads to a similar pile-up at very high densities, even in the MW & SMC cases (in which SNe heating alone can regulate global properties). 
Removing gas heating (by SNe and stellar winds) leads to more low-density cold gas and dramatically suppresses the hot gas mass. 



being smaller by a factor ~ 8, even though their instantaneous lu- 
minosity in the OB phase can be higher). 

The large variation in the amount and average density of the 
dense gas in Figure [12] and the deviations from log-normality, are 
striking, given that we showed in Figure [9] that the different feed- 
back models all produce nearly identical velocity dispersions and 
self-regulation at 2 ~ 1- They all maintain large-scale turbulence, 
in other words, even though some experience runaway collapse 
(which is then reflected in the much larger SFRs for these mod- 
els). But it is commonly assumed in various star-formation models 
that maintaining large-scale turbulence alone is sufficient to pre- 
vent runaway collapse (e.g.|Ballesteros-Paredes et al.|2007[|Tasker| 
|& Tan|[2009i [Krumholz & McKee||2005| l. Moreover many high- 
resolution studies of idealized, turbulent boxes have found that 
"pumping" the global modes in the box with momentum equivalent 
to 2 ~ 1 establishes a lognormal density PDF (with a dispersion 
that scales very weakly with mach number oc ^ln(l + 3M-/4); 
e.g.|Vazque z-SeiTia deni|1994[|Padoan et al.|1997||Scalo et al.|1998l 
lOstriker et al.|1999nFederrath et al.|200 8', 'Price et al."2011'l and 
prevents runaway collapse/star formation tV azquez-Semadeni et al. 
[20031 |Li et al.|2004||Li & Nakamura|20"06} . 

However, there are two major differences between our sim- 
ulations and many idealized turbulence models. First, the typical 
idealized box simulation realizes only a sub-region of one cloud: 
with a mass of a few tens typically measured in units of the thermal 
Jeans mass (~ G^^^^ p^^^^)\ as such the statistics only extend to 
the ~ 10% level (i.e. to ~ 1 a in the distribution), and "lognormal- 
ity" is seen only on a linear scale. Indeed, even in these simulations, 
hints have often appeared of non-lognormal "tails" in the distribu- 
tion (see ^Ballesteros-Paredes et al.|2011[ l. But with such sampling, 
it is not possible for a large sub-region (itself containing many Jeans 
masses) to detach from the global flow, dissipate, and collapse. In 
contrast our simulations here sample > 10^ thermal Jeans masses, 
and so can easily probe such fluctuations. 

Second, we include cooling and self-gravity, which are not 
always followed self-consistently in idealized simulations. This is 
critical for "runaway." Turbulence does tend to drive the distribu- 
tion to a lognormal. But the tails converge very slowly (with the 



Poisson error) - according to the central limit theorem, conver- 
gence to lognormality in a given portion of the distribution requires 
a number of "events" such that A'^ ' ^ 1 , i.e. for each Jeans mass 
region (here more properly defined in terms of the turbulent Jeans 
mass), A'c,.„ijsin= times ^ 1- if cooling is fast relative to the dy- 
namical time, a self-gravitating region collapses on a single cross- 
ing time. Outside of ~ 1 — 2 standard deviations from the core of 
the lognormal, then, a region cannot be "held up" or dissociated by 
cascades of turbulence from large scales "mixing" it back to low 
densities before it collapses. Regions will instead be continuously 
"pumped" to randomly cross this threshold, at which point they de- 
tach and collapse, building up the excess at large densities. It is 
therefore critical to invoke a feedback mechanism that can act di- 
rectly on the gas at very high densities to dissociate these collapsing 
regions before they completely run away, and re-mix them into the 
larger-scale medium. 

And this is exactly what we see: the fact that all the simu- 
lations maintain similar velocity dispersions is the reason that the 
median/peak density (around ~ 10 — lOOcm"^) in the dense gas 
is similar in all cases, and within ~ 1 cr of this peak, they almost 
all appear reasonably well-behaved lognormal distributions. But 
the radiation pressure and other feedback acting in the dense gas 
is critical to suppress the runaway "tails" seen in the no-feedback 
runs, and in the "No Radiative Momentum" runs in the three more 
massive galaxy models. It is these tails which in turn contribute the 
gas that is actually forming stars (densities > 10* cm~'), and hence 
lead to runaway star formation. 



8 THE PROPERTIES OF GMCS 

Having considered the global properties of the ISM above, we now 
turn to the sub-structure predicted: specifically, the properties of gi- 
ant molecular clouds. Figure [O] shows an illustrative example of a 
region near the center of the SMC simulation, with GMC structure 
resolved down to ~pc scales. It is clear that the GMCs formed ex- 
hibit considerable structure, with filamentary large-scale molecular 
gas and dense knots and clumps (the star-forming regions) embed- 
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Figure 13. Zoom-in of regions with a large concentration of GMCs in the 
SMC simulation. As in Figure[T] brightness encodes projected gas density 
(logaiithmically scaled with a ft: 4 dex stretch) and color encodes gas tem- 
perature, with blue being T < lOOOK molecular gas, pink ~ 10"* — 10^ K 
warm ionized gas, and yellow > 10^ K hot gas. The filamentary nature of 
clouds, with dense cores (bright white-blue regions) and diffuse molecular 
material (diffuse blue regions), and < lOpc-scale structure are evident. 



ded throughout. The irregular, clumpy, and triaxial morphology of 
these clouds is very similar to typical observed GMCs. They are not 
"small disks," which is the typical case when "GMCs" are formed 
in simulations with weak feedback (a point we will discuss in more 
detail below). 

To identify individual clouds (especially given their complex 
morphology), we follow a procedure motivated by observational 
studies, but taking advantage of the full simulation data. Specif- 
ically, we consider each simulation output snapshot and apply 
the sub-halo finder SUBFIND, which employs a friends-of-friends 
linking algorithm with an iterative unbinding procedure to robustly 
identify overdensities (for details and tests, see [Springel et al.| 
[2001 The routine is modified to run on the disk gas (instead of 
dark matter) and with softenings and linking lengths appropriately 
rescaled for the mean gas densities within the disk effective radius. 
Visual inspection confirms that this correctly identifies the obvi- 
ous clumps, and binding criteria (discussed below) also support the 
identifications. Changing the linking lengths and other numerical 
quantities has a weak effect on most properties, but we explicitly 
note where the GMC properties are sensitive to the algorithm. This 
tends to occur in the most massive clumps in each simulation, as 
these are often located at global critical points (e.g. along spiral 
arms or bars) and have their own sub-structure, so deciding whether 
or not to break them into sub-units is sometimes ambiguous. This 
same uncertainty, however, applies to observational studies as well 
- we discuss its possible effects where they may be important be- 
low, but these few cases do not affect our qualitative conclusions. 



8.1 GMC Mass Functions 

Figure [14] shows the resulting GMC mass function, for our differ- 
ent galaxy models with feedback enabled. The results converge to 
a quasi-steady state, despite the fact that (we will show) individ- 
ual clumps are short-lived, so the mass function reflects continuous 
dissociation and re-formation. The shape in all cases is an approx- 
imate power law with a steep (exponential or log-linear) cutoff at 
some high mass, so one could fit the usual truncated power-law or 
Schechter functions: 



A'(>M) =A^o 



exp(— M/Mo 



(10) 



we choose this form so that a corresponds to the logarithmic slope 
of dN /dM, standard in observational studies. Observed values tend 
to lie within the range —2.1 ^ ol < —1.5 shown in the Figure [T4] 
with the "canonical" value of —1.8 very close to the typical slope 
we find (fitting the function above to all simulation timesteps, we 
find an average a ~ — 1 .75 ± 0.2). 

We show this explicitly by comparing the observed GMC 
mass functions for the MW (from [Williams & McKee|1997| cor- 
rected to represent the total MW population, since we consider all 
GMCs in the simulation), appropriate to compare to our MW-like 
simulations, and for the LMC (from [Rosolowsky [2005 | l, appropri- 
ate to compare to the SMC-like simulation. Of course the normal- 
ization A'o of the observed and simulated populations is arbitrary 
and depends on the total gas mass, but for these comparisons the 
total gas mass is similar. We expect the Sbc case to be qualitatively 
similar to these but with different normalization and zero point. For 
the high-redshift case, the total gas mass and most massive clump 
mass are much larger than local systems, so there is no direct com- 
parison available, but observations have indicated that typical sys- 
tems similar to our HiZ model have approximately ~ 5 — 10 clumps 
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Figure 14. Predicted mass function of GMCs in the simulations (see ^for a discussion of how we identify GMCs). We show the different galaxy models, 
each with runs that include different feedback mechanisms (see Table 2). The distributions can all be fit by a power-law with a cut-off at high mass (eqn.flOl. 
For the MW model and dwaif/SMC model, we compare with the observed MW and LMC mass functions, respectively ^Williams & McKee 1997 , Rosolowsky| 
|2005| orange diamonds). Characteristic slopes in obsen'ed systems are also shown for comparison (orange lines for a = — 1.5, — 1.8, —2.1, where A'(> M) oc 
M"""*"'). In all models, the predicted slope is similar to that observed: it follows generically from gravitational collapse (with super-sonic flows). The cutoff 
mass traces the Jeans mass in each model (eqn. |ll( ; since all reach Q ~ 1, the high mass cutoff is weakly sensitive to feedback, primarily through more/less 
efficient gas exhaustion (/gas{S)). With weak/absent feedback, the very high-mass tail is enhanced even though typical clouds are less massive; this is because 
massive clouds along global ISM structures (e.g., spiral waves) do not disrupt but accrete continuously. 



in the ^ l(f Mq range, similar to the number predicted (Forster 
|Schreiberetal.|2011t . 

The mass function slope appears fairly generic. To see this, 
we compare models with different feedback mechanisms enabled 
and disabled, as in §|6] The characteristic or maximum masses may 
shift by small amounts without certain feedback included, but the 
slope Q is similar in all cases. We have also checked this within 
each feedback model, for example, arbitrarily multiplying or divid- 
ing the strength of the momentum coupling from radiation pressure 
by a factor of 5. A mass function slope close to a ~ —2 tends to 
emerge generically in any system dominated by gravitational col- 
lapse (e.g. [Press & Schechter|[T974i [Bond et al.|[T99T] l; the loca- 
tion of the high-mass cutoff and exact deviation from —2 depend 
on details of e.g. the power spectrum, non-gravitational terms, and 
Jeans conditions (if the medium is gaseous), but the low-mass slope 
must be close to —2 because gravity is scale-free[^This explains 



° Implicit here is the assumption that motions are super-sonic, so the ther- 
mal properties do not set a prefeired length scale, but this is easily satisfied 
in the simulations and observations, down to scales well below what we 
resolve (< 0.03 pc). 



why simulations with very different physics included recover simi- 
lar behavior (compare our runs and e.g. [Audit & Hermebelle|20T0} 
|Tasker|201l) . 

The cutoff mass Mq appears to be simply related to the turbu- 
lent Jeans mass in each disk: 



Mj = =4TvT,h^ 



(11) 



where the latter expression includes the Toomre 2 ~ 1 and the pa- 
rameter v—l — l which depends on the mass profile, and /gas refers 
to the fraction of gas mass relative to the total enclosed mass. Eval- 
uated at Rg for the typical conditions of each galaxy model seen 
in Figure [9] this gives ~ (1 - 5 x 10*, 0.5 - 2 x lO'', 1 - 4 x 
10', 0.3 - 1 X 10'')Mq for the SMC, MW, She, and HiZ cases. Of 
course, the exact values will change with time and the precise loca- 
tion in the disk, and mergers and agglomerations can grow clumps 
beyond this limit. But it appears to be a quite good approximation 
to the "cutoff" in the integrated mass functions. 

When we compare models with different feedback mecha- 
nisms enabled or disabled, the slope is similar, but there are some 
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small shifts in the "cutoff" mass. At low/intermediate masses, mod- 
els with less efficient feedback are shifted to slightly lower GMC 
masses. This is most noticeable when we remove radiation pres- 
sure in the HiZ or She cases, remove hot gas in the SMC case, or 
remove all feedback in any case - just what we expect, since we 
already saw these have the largest effects in Figure [9] Recall, we 
showed above that all the models tend to 2 ~ 1, so this primarily 
reflects somewhat more efficient gas exhaustion lowering Mgas and 
so the Jeans mass. 

However at the very highest masses the behavior is perhaps 
unexpected; despite these shifts the extreme high-mass "tail" of 
the mass function extends to higher maximum masses. Two effects 
drive this. First, the most massive clumps tend to be massive clump- 
complexes along e.g. spiral arms, which are linked by intermediate- 
density material. Without being able to fully dissociate this mate- 
rial, the friends-of-friends GMC finder tends to link these into sin- 
gle "super-clumps." This alone is ambiguous - it is well known in 
observations that at high masses the assignment of "clump mass" 
is quite sensitive to the exact clump identification and linking cri- 
terion (see [Pineda et al.|2009[ l. We err on the side of over-linking 
rather than missing real clumps, but if this alone drove the effect 
seen we would not assign it much meaning. However there is also a 
robust physical effect: clumps survive longer with weak feedback. 
The same most massive clumps, located along global structures, 
can then accrete large quantities of gas, growing beyond the Jeans 
mass. This leads to the odd shapes of the high-mass mass function 
in these cases, where the high-mass end does not rapidly fall off, 
but "turns up" in a manner not observed to date. 



8.2 Density Distributions 

8.2.1 Three-Dimensional Densities 

The distribution of three-dimensional densities n in clouds is 
more or less the distribution of "cold gas" shown in Figure [T2] 
As discussed above, the distributions in our standard models are 
crudely log-normal with median ~ 10 — lOOcm^^ and scatter 
0.25 — 0.30dex, within which most of the GMC mass resides. 
This value increases systematically with the average density of the 
galaxy; fitting a mode to each lognormal gives « 25cm~' in the 
SMC model and lOOcm"^ in the HiZ model. This must happen: 
the HiZ model has an initial smooth disk-averaged n ~ 40cm~^ 
inside Re, so any clumps must have larger n. There are tails at the 
< 10% level in the distributions which include the high-n cores. 
This reproduces a key feature of observed GMCs: although most 
of the mass within e.g. the star-forming galactic disk is in "dense" 
gas, in the form of GMCs, most of the gas within those GMCs is 
at relatively low densities, not within star-forming cores that have 
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densities > 10 cm~ (e.g. 
and references therein). 

As noted above, these distributions remain steady-state 
throughout the simulations. After the initial couple dynamical 
times, the median densities vary (randomly) within factors < 1.5 
and dispersions vary within 0. 1 dex. This is non-trivial - as time 
passes, there is not a one-way collapse to high densities but self- 
regulation at a constant density distribution. 

However, models without feedback do not have a stable den- 
sity distribution, and "pile up" at the highest resolvable densities 
n S> lO^cm"' as seen in Figure 
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We saw there that this is pri- 
marily a function of the radiation pressure support. In the MW and 
SMC models, gas heating by stellar winds and HII photo-ionization 
is able to stave off complete runaway cloud collapse but without ra- 
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Figure 15. Distribution of average GMC surface densities (S), calculated 
as the mass-weighted (E) averaged over viewing angles for each cloud. We 
show each disk model with different feedback mechanisms enabled (see 
Table 2). The distributions are close to log-normal in the core, with the 
median surface density ~ 50 — 300 Mq pc~^ increasingly modestly as the 
average gas surface density of the star-forming disk increases. The scatter 
in the surface density of GMCs is small (^ 0.2 dex). For comparison, we 
show the observed median (S) and it 1 tr range (shaded orange) inferred for 
SMC and MW clouds and the range for all dwarf' galaxies observed for the 
Sbc-like run (Bolatto et al. 2008 |; for the HiZ run we take the estimate for 
4 ~ 2 galaxies in Forster Schreiber et al. If 2009 1. With no feedback, local 
collapse proceeds without limit (to high n), which gives rise to the high-S 
tail in the GMC surface densities seen here, but does not necessarily change 
the median (S) . With feedback, clouds are more triaxial and short-lived and 
so do not globally contract by a large amount. 



diation pressure the number of ultra-dense cores is enhanced. We 
emphasize though, that the effect of feedback is primarily on the 
high-density tails in Figure[T2j since most of the gas is at low (non 
star-forming densities), its distribution is similar even in the pres- 
ence of runaway local collapse and star formation. 



8.2.2 GMC Surface Densities & Size-Mass Relation 

A closely related metric is the average cloud surface density (S), 
shown in Figure [TS] Here we use the mass-weighted average sur- 
face density over sightlines through the cloud (E) ; this is not ex- 
actly the same as the area-average ~ M/tt^?^ but is well-defined 
even for arbitrary cloud geometry. In each model, the (E) distri- 
bution follows a relatively narrow log-normal. The median GMC 
surface density increases from ~ 50Mqpc~^ for the SMC, to 
~ 500 Mq pc^^ in the HiZ model, as they must given the increas- 
ing background density. The scatter is again roughly constant, only 
~ 0.2 dex, and the distributions in this projection are closer to exact 
log-normal than n. This narrow distribution in E means the clouds 
follow a size-mass relation approximately R oc M""^, with a normal- 
ization consistent with the observed range (weakly increasing with 
galaxy surface density) and small scatter. 

These surface densities (and corresponding size-mass rela- 
tion) and their scatter are quite similar to those observed in MW- 
like and SMC-like local galaxies ([Bolatto et al.|2008"i|Fukui et ah] 



[2008i|Wong et al.|2008[[Goodman et aL|2009| l. Observations have 
even found evidence for a systematic shift in the normalization of 
the size-mass relation ((E)) towards lower densities in SMC clouds 
and higher densities in dense regions of massive galaxies jBolatto[ 
[et al.|2008| [Heyer et al.|2009| l. In fact, the properties of the most 
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massive clouds in our HiZ simulation agree well with those inferred 
for the most massive star-forming molecular "clumps" in z ~ 2 — 3 
galaxies observed with integral field data and in some cases individ- 
ually lensed jForster Schreiber et al.|2009{|Swinbank et al.|201 
which are offset towards higher (E) relative to the extrapolation of 
the local relation. 

Our predictions are also, interestingly, similar to the predic- 
tions in higher-resolution simulations of individual GMCs with 
quite different physics included (e.g. ignoring radiation pressure 
but including magnetic fields; see [Ostriker et al.|[2001[ |Tasker| 
|201 1\ . This suggests that they are fairly generic. No-feedback cases 
do feature a larger tail towards high-E, but the median values are 
actually quite similar to the cases with feedback. This fits with the 
behavior we saw in § |6.1| where dissipation is efficient in the ver- 
tical direction of the disk without feedback, but there is still large 
radial and azimuthal motion. In detail, what we see with no feed- 
back is often clouds collapsing, dissipating all of their energy in 
a crossing time, but then being arrested by residual angular mo- 
mentum. They become small, rapidly rotating disks, with clumpy 
sub-regions collapsing further. This has been seen in other high- 
resolutions disk simulations without feedback, as well l |Boumau^ 
[erari[2007l [Pobbs 2008, [Tasker & Tan 2009). In other words, 
without feedback, GMCs do not collapse three-dimensionally, and 
so they can locally collapse to arbitrarily high densities, turning 
all their mass into stars in a couple crossing times, even though 
their surface densities are not very different from cases with ef- 
ficient feedback. This behavior, and strong rotational support for 
highly "flattened" molecular clouds are, however, not observed 
( [Rosolowsky et af]|2003[ |Imara et al.||201ll l. Much closer to ob- 
served systems are the models with strong feedback, which have 
similar local and aR, leading to molecular clouds which, while 
certainly non-spherical, are not highly flattened (they tend towards 
irregular, triaxial and filamentary shapes, with some preference for 
elongation in the galactic plane as observed; see jKoda et al.|2006) 
and not rotationally supported (see Figure [T3]l. 

With feedback included, clouds are three-dimensional; feed- 
back has greatly slowed collapse. Moreover, we will show that with 
feedback, clouds are dissociated in just a few dynamical times. As a 
result, the clouds can undergo little global contraction (even though 
over-dense regions within the clouds can collapse to very high den- 
sities and form stars). The surface densities (E) are therefore just 
a reflection of the average surface densities of the disk from which 
they form. Specifically, the peak in typical cloud surface density 
in Figure^] in each model is just ~ 2 — 5 times the average disk 
surface density within the star-forming disk. In other words, the 
remarkable independence of the global cloud densities on the de- 
tails of feedback, and apparent narrow range in densities, is just a 
consequence of the fact that they are dissociated after just a few 
dynamical times, before they can globally contract (and are slowed 
in that contraction even during that time). This is quite similar to 
what is observed: although the star-forming cores within GMCs 
reach large columns Nb > 10^''cm~^ (optically thick in the IR), 
MW GMC complexes are observed to have a surface-averaged col- 
umn density of just 10^' — 10^^ cm~^, similar to the average through 
the disk. 

8.3 Linewidth-Size Relation 

Closely related to the typical densities and surface densities of 
clouds are the Larson's law-type scaling relations. We have al- 
ready discussed the size-mass relation; in Figure [16] we plot the 
line width- size relation. For simplicity, we define the radii and ve- 



locity as the rms spherical radii and velocity dispersion; because 
observations generally de-project to three-dimensional quantities 
assuming an ellipsoid in projection, this is probably a reasonable 
approximation at the level we consider here. We are not making 
a detailed mock observation of molecular gas in various states, so 
there will be some systematic uncertainties in the comparison, but 
the senses of the scalings are robust to various density limits and 
definitions. For clarity, we plot just a small randomly sampled sub- 
set of the clouds. 

The predicted linewidth-size relation is broadly similar to that 
observed with an approximate a oc I^ ^ scaling, in the range where 
there is overlap (~ 5 — 300pc), and continues the observed scal- 
ings beyond that to larger ~kpc clump sizes in the massive sys- 
tems (we plot two characteristic points from the observations of 
IForster Schreiber et aI.|2009||F6rster Schreiber et al.|201 1| for z ~ 2 
clumpy disks). As expected from the behavior of (E), there is a nor- 
malization shift such that the HiZ systems obey a similar power-law 
but with higher a at fixed R (reflecting their higher densities). The 
predicted values, which exhibit this shift, appear consistent with 
the observed massive high-redshift systems. For clarity, we have 
shown just our "standard" model, because the power-law scaling 
is similar for most of our models, with normalization offsets that 
follow simply from the differences in average density (above) and 
virial parameters (discussed below). 

Again, models with a wide range of different physics included 
appear to capture the linewidth-size relations. As with the size- 
mass relation, the key point is that with feedback present, clouds 
are short-lived and collapse weakly. Moreover because the scaling 
above is based on second moments in velocity and size, it is es- 
pecially sensitive to the less tightly-bound/collapsing material. As 
such, the observed scaling reflects the conditions of marginal self- 
gravity and Jeans collapse. The scale-length of Jeans collapse is 
just A ~ (Sv^/ttGE. So if the systems have not been able to evolve 
very far from their initial collapse, i.e. if they are short lived, then 
we simply expect Sv oc \/A across regions with similar E , with a 
residual normalization dependence 



similar to what is observed and simulated. 
8.4 GMC Virial Parameters 

Figure[T7]shows the distribution of virial parameters a of the simu- 
lated clouds, defined as the ratio of kinetic to potential energy. Ob- 
servationally, this must be estimated from global cloud properties, 
so is commonly defined as 

" ~ I V I ~ fl2 GMl/R,x ^ GMci Mci 

where the factor of 5 follows from common assumptions for the 
terms ai and 02 (order-unity constants that depend on the true mass 
profile shape and projection), a and Rg are the average dispersion 
and radiusrland Mvi, = 5a~ R/G. With this definition (if the as- 
sumed coefficients are correct), a = 1 would be virial equilibrium. 

The "average" clumps are consistent with being marginally 
gravitationally bound, a ~ 1. There is a broad dispersion of > 
0.5 dex about this median, reflecting both different states of various 

' For consistency with the observations we compare to, we take a and Rg 
as randomly-sampled line-of-sight and corresponding projected dispersion 
and half-mass-radii. 
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Figure 16. Linewidth-size relation of GMCs in our standard simulations of each galaxy model. We randomly sample clouds in time for each disk model. We 
only show clouds with > 100 member pa rticles and sizes > 5 softening lengths (this determines the low-mass cutoff in the plots). We compare observations 
from local-group galaxies [Bolatto et aE H2OO8. filled points), with typical best-fit power-laws a oc R"'^ (dotted lines). For the HiZ case, we compare typical 
properties of very massive clumps observed in |Forster Schreiber et al.|j2009^ ; |F6rster Schreiber et al.|^20lT) (open points). Because GMCs are short-lived, 
these scalings essentially reflect the initial collapse conditions, and so are relatively insensitive to feedback. 




Log[ Virial a = 2T/|V| ] 

Figure 17. Predicted distribution of GMC virial parameters (the ratio of 
kinetic to potential energy; eqn. |13) . We plot the distribution for all disk 
models together, since the individual results are very similar Once sufficient 
feedback exists to resist collapse and make clouds relatively short-lived, all 
models equilibrate at marginal binding (a ~ 1 and is only weakly correlated 
with feedback strength); With no feedback, however, clouds are much more 
tightly bound, with global a < 0.1. 



clouds, but also projection effects, since the clouds are highly non- 
circular. This is all expected, since we argued above that clumps 
dissociate before they can become strongly bound or globally con- 
tract much. 



Within each cloud, there is a broad dispersion as well of "lo- 
cal" virial parameters. Much of the diffuse material at (n) is actu- 
ally un-bound (at least ~half the material is at a > 1, and can reach 
a ~ 10 — 100). The bulk velocity of this material, in the center of 
mass frame of the GMC, exceeds the cloud escape velocity. This 
material is associated with the clump only in a transient sense (e.g. 
from turbulent compression, shocks, outflows from dense regions, 
or convergent flows). But dense clumps with local n > 10''cm~^ 
(the regions that will actually form stars, but a small mass fraction) 
tend to be bound. This is similar to what has been inferred from 
observations of GMC complexes (e.g. |Heyer et al.|200 1 1. Quanti- 
tatively, we can consider the "true" Otme (calculated knowing the 
actual true kinetic and potential energy relative to the cloud cen- 
ter of mass for every gas particle) for two separate populations: 
the cloud "cores" (for convenience, gas inside the spherical half- 
mass radius with relative velocity below the median of the clump 
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Figure 18. Correlation of the GMC virial parameter a (Figure [TT) with 
mass for our standard simulations with all feedback mechanisms enabled. 
We compare the observed local group clouds from |Solomon et ar]jl987 |. 
The model and observations both show a weak trend in which more massive 
GMCs tend to be more tightly bound a oc M^"^. 



material), and the remaining "non-core" material; we find typical 
(atrue) ~ 0.05-0.2 and («,rue) ~ 30- 100, respectively!^ 

With no feedback, clouds are at least an order-of-magnitude 
more tightly bound, a --^ 0.05 — 0. 1 . This follows from their (largely 
vertical) collapse, giving smaller R and a at the same mass (these 
quantities often collapse together in a way that keeps the clouds 
not far off the linewidth-size relation, but gives much smaller a 
at fixed mass). However, once sufficiently strong feedback is in 
place to slow/dissociate clumps, the properties are relatively insen- 
sitive to the feedback model. This is not surprising: provided suffi- 
cient energy/momentum to slow Jeans collapse and/or dissociate a 
cloud, objects will either globally equilibrate at marginal binding, 
or spend most of their lifetime at the point where infall reverses to 
expansion, i.e., at the maximum value of a. If the feedback operates 
rapidly, this minimum a will be near the value where star forma- 
tion commences. This is why (as we saw with other cloud proper- 
ties) models with very different assumptions and more simplified 
models for feedback have also seen similar characteristic a ( |Dobbs| 
|et al.|20l"T| l. For the HiZ case, we again find that radiation pressure 
momentum is most important (removing sources of hot gas has lit- 
tle effect). For the SMC-like case, in contrast, gas heating plays a 
large role in regulating the virial parameters. 

Figure[T8]plots the virial parameters of GMCs as a function of 
their mass. We show this only for our standard model, as the sys- 
tematic offsets between models can easily be rea d off Figure | 17| 
For comparison, we plot the observed points from [Solomon et aT] 
( |1987^ . The simulated and observed GMCs agree well. Not only is 
the median a similar, but the simulations reproduce the observed 
weak scaling of a with Md, approximately a oc Mj"'', such that 
more massive clumps tend to be more bound (see |Heyer et al.| 
|2001[ l. This is partly related to the fact we show below, that more 
massive clumps tend to have longer lifetimes. We note, however, 
that Figure[T8]implies that large fractions of the population are re- 
ally not self-gravitating - it is better in this regime to consider the 



systems as simply "molecular overdensities" rather than "clouds" 
in the traditional sense. 



8.5 Lifetimes and Star Formation Efficiencies 

We next consider typical cloud "lifetimes" and integrated star for- 
mation efficiencies. To define these, we need to link clouds in time 
between different snapshots. Given a specific clump p in snapshot 
i, we identify the descendant d of this clump (in snapshot ; + 1) 
as the clump which contains the most total mass in particles from 
the original clump p. If no clump in snapshot / + 1 contains any 
particles from p, then the clump has no descendant (its mass be- 
comes zero). This defines a "merger tree" for clumps with timep^ 
To consider a "main branch" of that tree (i.e. growth of the pri- 
mary clump, which may grow by accretion/inflow, or by mergers 
of smaller clumps), we note that if a clump d is identified with 
multiple progenitor clumps p, the "main progenitor" is simply the 
most massive. The other progenitors are marked as having "merged 
onto" the main branch at this time. 

For each main branch, we then have a mass evolution versus 
time, Mci(f). This generally rises exponentially as the clump starts 
forming, then (with feedback enabled) falls rapidly as feedback un- 
binds the gas. We can then define other important quantities: the 
maximum mass of the clump (just M^iax = MAX(Mci[r])), and the 
lifetime fci- We take the latter to be the total time when the clump 
is above some cutoff threshold rj ~ 10% of M,„ax - the choice is ar- 
bitrary, but because the clumps tend to grow and dissociate quickly 
the lifetimes are not especially sensitive to the choice so long as 
??< 1. 

Figure [19] plots the distribution of the resulting clump life- 
times in our standard simulations. We plot it in absolute units and 
in units of the clump dynamical time, fdyn = 1 / \/Gpci- We do not 
include clumps that are lost via merger onto more massive clumps, 
because we cannot know how long they would have survived, but 
their growth/decay curves tend to be no different from main-branch 
clumps of similar mass and spatial locations until their merger. And 
we find that the majority of clouds grow primarily by accreting 
"smooth" gas (i.e. gas not in another massive cloud), rather than 
via (at least major) hierarchical mergers, so this is not a large ef- 
fect. 

Our simulated GMCs are short-lived: with feedback present, 
most live ~ 10* — lO'yr, ~ 1 — 10 free-fall times. Star formation 
within clouds should therefore be inefficient: we expect just a few 
percent of their mass will be turned into stars, but we can calcu- 
late this explicitly. We sum the SFR from all cloud particles at 
each time, integrated over the cloud lifetime to get the total stellar 
mass formed, and compare this to the maximum clump mass (defin- 
ing the efficiency M,.fonned/A^cioud,max). The typical cloud converts 
~ I — 10% of its peak mass into stars, similar to what has been 
inferred from a range of observations i Zuckerman & Evans|1974[ 



[Williams & McKee||l997l [Evans[|T999[ [Evans et al.[[2009| >. This 
does not necessarily mean, however, that this fraction of all mass 
that enters the cloud is converted into stars. Over the cloud lifetime, 
mass is continuously accreted and lost, so the total mass in gas that 
gets processed "through" the cloud can be even larger, implying an 
even lower net "efficiency" of converting GMC material to stars. 



Note that because the "diffuse" material is a non-negligible fraction of 
the mass, calculating the "true" a for all cloud material often gives large 
values, dominated by the material which is transiently associated in shocks 
and has large kinetic energy. 



' ' To properly sample this requires that the snapshot time spacing be less 
than typical clump lifetimes. We have experimented with spacings as short 
as lO' yr, and generally find converged results for the lifetime statistics pre- 
sented here for spacings < 10^ yr. 
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Figure 19. Statistics of cloud lifetimes and star formation efficiencies in our standard simulations with all feedback mechanisms enabled. For each simulation, 
we identify all clumps that can be well-resolved (> 100 particles) down to 10% of the clump mass (> 1000 particles at peak); lifetime is defined as the time 
above 10% the maximum cloud mass. Top Left: Distribution of cloud lifetimes for each of our galaxy models. Top Center: Lifetimes relative to the free-fall 
time t^x/ {t[f) of each cloud. Top Right: Integrated star formation efficiency (stellar mass formed in clump over maximum cloud mass). Bottom Left: Lifetime 
versus maximum cloud mass. Bottom Center: Stellar mass formed versus maximum mass. Lines correspond to the different efficiencies labeled. Clouds live 
for a few dynamical times, and turn a few percent of their mass into stars. There is also a weak trend for more massive systems to have higher efficiencies and 
longer lifetimes. Bottom Right: Integrated star formation efficiency (as top right), but in models with no feedback. Without feedback, clouds persist until most 
of their dense gas is turned into stars, giving star formation efficiencies of ~ 10 — 50% on timescales of ~ lOMyr. 



What sets the efficiencies and lifetimes of GMCs? As we dis- 
cussed above, for the dense gas, the most important feedback mech- 
anism is the local radiation pressure. Therefore consider just this, 
for simplicity, in a spherically symmetric cloud of mass M and pro- 
jected effective radius Re. The total gravitational force in a GMC 
(summing over mass) is 



GM- 



(14) 



where /3i ~ 1 depends on the mass profile The force from radiation 
pressure, integrated over the sphere, is 

L 



i^rad = (1 +-riR) 



(15) 



where the pre-factor assumes the cloud is optically thick in the 
UV/optical, but not necessarily in the IR. If the lifetimes are short, 
of order a few Myr, then L oc M, , because most SNe have not ex- 
ploded; for a standard stellar population with the same IMF as- 
sumptions as our models, L ~ 1 137 Lq {M,/Mq). Since the aver- 
age surface density of clouds is within a narrow range, it is con- 
venient to define (E) = I32M/{ttRI) ((E) is the mass-weighted 
average E, ^2 again just depends on the mass profile). Setting 
^'gmv = ^'rad, We obtain 



Mclo, 



!0.04/3o.7 



1+riR 



(16) 



/0.7 (= 0.95, 1. 07, 1. 12 for a constant- 
Hemquist 1990 profile) and {E)ioo ~ 



where /3o.7 = (Pi ^2)- 
density, Plummer, or 
(E)/100Mqpc--. 

This simple argument reproduces the efficiencies at low cloud 
masses, where these assumptions hold, and tir is small, for the typ- 
ical observed average surface densities of such clouds ((E) 100 ~ 1). 
The predicted efficiency does increase at higher (E), which could 
in principle lead to a runaway collapse if there were no offset- 
ting force. However, at sufficiently high E > IOOOMqpc"^, tir 
becomes larger than unity even for the smoothed average surface 
density. More accurately the typical tir ~ c" Km (E) where c~ is a 
"clumping fudge factor" that accounts for the excess rise in density 
around each stellar cluster/core sub-region, a well as contributions 
from other feedback mechanisms (typical effective c ~ 2 — 5). 

Once riR ^ 1 , the efficiency should asymptote to a maximum 



■0.4 c -pojKf 



(riR > 1) 



(17) 



where ^5 — KiR/5cm^g~' ~ 1. Given modest c ~ a couple, and 
the fact that our simulations rarely reach such extreme densities, 
this matches reasonably well the "upper limit" to the star formation 
efficiencies we find; together, these values bracket the range in the 
simulations. 

Cloud lifetimes follow from these integrated efficiencies, as 
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the time needed to form this number of stars, Af = M, / M, (cloud) . 
For an observed instantaneous SFR of a one or two percent per free- 
fall time, it takes a few free-fall times to form this mass and disrupt 
the cloud (in this limit, of course, we insert this efficiency in our 
SFR model, so the lifetime is not predicted in a fully a priori sense). 
However, in Paper I, we showed that the global star formation rate 
did not depend on the efficiency employed in our model. Since it 
takes a free-fall time to assemble a GMC, and a free-fall time to 
disrupt it, even if we assume a large efficiency, the cloud will live 
one to two free fall times. 

In the most massive cloud complexes, the lifetime reaches 
as large as > 30Myr and is no longer short compared to stel- 
lar evolution timescales. Since massive stars are dying, instead 
of L oc M,, we have L oc M, . In this limit, for a constant star 
formation rate L/M, ^ 5.7 x lO' Lq/(Mq yr~'). Inserting this 
into the derivation above requires an instantaneous SFR M, ~ 
CSMgyr-' (Mcm/ lO' Mq ) /3o.7 {E),oo (1 + tir)"'. This then re- 
CJUirCS 3. tiniG ^lifetime ~ fdyn to dissociate the cloud (and must act 
in less than a few dynamical times so the momentum is not dissi- 
pated); if we require the total momentum injected at this rate equal 
the binding momentum of the cloud, we obtain fufctimc ~ 1.4/3'rff 
(where again depends very weakly on the mass profile and 
(te) = (37r/32{p))'/^), giving 



^« 0.04/3' " 

Mcloud 1 + TIR 



4Myr/ 



;0.04^" 



1+riR VlOpc/ 



(18) 
(19) 



Note that for similar-{E) clouds, the characteristic dynamical 
time increases with cloud mass as fff oc ^J? oc M'^j^^^. We there- 
fore expect more massive clouds to have larger absolute lifetimes, 
a correlation clearly evident in Figure [T9] Moreover, according to 
Equation [TS] the integrated star formation efficiency (in massive 
clouds, at least) should also be a weakly increasing function of 
the cloud mass oc mJ?io„j. We also see this trend in our simulations, 
with efficiencies rising from 1% to 10% from M,, < IO^Mq 

to Mci > 1 0*^ Mq , or M« . formed /M doud. max OC M^°^l^°^2 (^se also the 
discussion i n | Murray et al.|2010| l. 

Figure [19| compares the cloud star formation efficiencies in 
runs without feedback. In absolute units, the lifetimes are not nec- 
essarily much longer - many live ~ lOMyr; however the popula- 
tion with lifetimes ~ 1 Myr disappears and a large tail of clouds 
(often in global structures such as spiral arms) appears with life- 
times ~ lOOMyr. But, recall that in the feedback-free cases, the 
dense gas collapses rapidly. As such, the range of lifetimes in units 
of the free-fall time oc p^'''' is much more broad, with clumps liv- 
ing as long as --^ 10 — lOOfe. The distributions are very broad be- 
cause it is often local regions within the clump running away to 
high densities, so the average (p) may or may not reflect this to 
varying degrees in a cloud. The more robust indicator is the inte- 
grated star formation efficiencies. As expected, from the runaway 
global SFR in these models, the star formation efficiency per-cloud 
is ~ 10 - 100%, with a mean of ~ 40% (it is not 100% because 
there is always some low-density material associated with clouds, 
which is left at even lower densities once dense regions have all col- 
lapsed). Feedback is necessary to reproduce the observed per-cloud 
cloud lifetimes and star formation efficiencies. 

Considering each feedback mechanism in turn, our results es- 
sentially mirror those for the cold gas and GMC density distribu- 
tions. In all models, radiation pressure is important to resist col- 
lapse; removing it in the HiZ case gives results nearly identical 



to the feedback-free case. And runs in all models with only radi- 
ation pressure enabled give similar results. In the MW and SMC 
cases, gas heating on ~ 1 Myr timescales from "fast" stellar winds 
and photo-ionization can contribute comparably to cloud dissocia- 
tion (removing radiation pressure does not result in total collapse if 
they remain, though the average star formation efficiency is slightly 
higher), and the longer-lived clouds in these models tend to be dis- 
sociated when SNe turn on (a few Myr). 

9 DISCUSSION 

We have implemented a model for stellar feedback in SPH simula- 
tions of galaxies that incorporates physically motivated treatments 
of a number of the key feedback processes: radiation pressure from 
UV, optical, and IR photons; the momentum injection, gas heating, 
and gas recycling from SNe (Types I & II) and stellar winds; and 
HII photoionization-heating. Our treatment of radiation pressure is 
calibrated by comparison to more detailed radiative transfer calcu- 
lations (Appendix [A]i. Our model of the ISM also includes molec- 
ular cooling, a simple treatment of HII destruction, and allows star 
formation only at very high densities within giant molecular clouds 
(n > 100cm~^^). 

In this paper, we study the consequences of these feedback 
models for the global structure of the ISM in disk galaxies (veloc- 
ity dispersions, Toomre Q, etc), the phase structure of the ISM, and 
the properties of GMCs. We simulate four representative galaxy 
models to bracket the diversity of star-forming galaxies observed 
locally and at high redshift (Table[2j Figs. |l|5| ). In a companion pa- 
per, we study the galaxy-scale outflows driven by stellar feedback 
and show that galactic winds appear to remove the bulk of the gas 
from their host galaxies. 

With physically plausible stellar feedback mechanisms en- 
abled, we show that star-forming galaxies generically approach 
a quasi steady-state in which the properties of the ISM are self- 
regulated by the turbulence driven by stellar feedback. GMCs form 
by self-gravity and survive for a few dynamical times before they 
are disrupted by stellar feedback; during this time, GMCs turn a 
few percent of their mass into stars (Fig.|19|>. The resulting galaxy- 
integrated star formation rates are in reasonable agreement with 
the Kennicutt-Schmidt law (Fig.|7j see also Paper I) and the global 
galaxy properties (dispersion, scale heights, Toomre Q parameter; 
Fig. [9j and phase distribution (Fig. [TTJ of the ISM are consistent 
with observations. 

The ISM phase distributions calculated in this paper are of 
course still limited in accuracy both by standard numerical restric- 
tions and by the inability of any galaxy-scale calculation to capture 
small-scale mixing processes and the effects of (saturated) thermal 
conduction. Nonetheless, much of the ISM dynamics and phase 
structure we find is a consequence of supersonic turbulence driven 
by stellar feedback. Such turbulence is reasonably well modeled 
numerically by SPH (based on comparisons to grid-based meth- 
ods; e.g.,|Kitsionas et al.|200^ [Price & Federrath|2010[ [Bauer &| 



[Springer 201 1). Based on this and preliminary comparison with 
other codes, we expect that the inclusion or exclusion of different 
feedback mechanisms will result in much larger differences than 
the details of the numerical method. 



9.1 The Problems of Feedback-Free Models 

The structure of the ISM in the presence of stellar feedback is 
markedly different than in calculations that do not have any stellar 
feedback. In the latter case, the gas rapidly radiates away its vertical 
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thermal support and any gravitationally-induced turbulent fluctua- 
tions are damped on a crossing time. This leads to runaway col- 
lapse of the gas and star formation well in excess of the Kennicutt- 
Schmidt relation (M* ~ Mgas/fdyn; Fig.^. Dense star-forming re- 
gions within GMCs collapse without limit, producing a large non- 
lognormal excess of high-density gas (Fig.|12|(. These GMCs have 
very low virial parameters relative to observations (Fig. and 
turn ~ 10 — 50% of their gas into stars into just a few cloud dy- 
namical times (Fig.|19|l. 

In other words, turbulence from global gravitational instabili- 
ties and the collapse of GMCs does not appear sufficient to main- 
tain low star formation efficiencies and prevent runaway collapse of 
the ISM (neither globally nor within GMCs). This result has been 
found in other simulations with sufficiently high resolution to fol- 
low individual GMCs (see e.g. |Boumaud et al.|20l"0l|Tasker|20TT| 
|Dobbs et al.|20TT] l. Although our calculations do not include inflow 
of gas from the IGM and thus we cannot explicitly assess whether 
such accretion can help maintain turbulence in galactic disks, the 
fact that the runaway is local, and occurs even while the disk main- 
tains a large nominal gas velocity dispersion, suggests that gas in- 
flow will not change these conclusions. 

We caution that some properties of the ISM can be reasonable 
in feedback-free models, but for the wrong reasons. For example, 
without feedback disks still eventually reach (2 ~ 1 (Fig.|9](; this 
is true even without star formation so long as cooling is included. 
However, without feedback this occurs because of two processes 
that do not occur in real galaxies: first, all of the dense gas rapidly 
turns into stars, leaving only a small, tenuous gas mass behind that 
can more easily maintain 2 ~ 1 ; second, in doing so, the dense gas 
collapses without limit to > 10* cm^' regions that become so small 
that they dynamically act like the coUisionless stellar disks. 

In our simulations without stellar feedback, the global col- 
lapse of GMCs is eventually arrested by angular momentum f*^ The 
resulting GMCs have roughly similar median surface densities and 
Larson's law scalings to models with feedback (Fig. [Tsj, but dy- 
namically the GMCs are rotationally supported disks instead of tri- 
axial clouds. Moreover three-dimensional collapse still runs away, 
especially in dense sub-regions. Thus the median surface densi- 
ties of GMCs are similar with and without feedback, but the mass 
in dense molecular gas and the per-cloud SFRs are an order-of- 
magnitude larger in the absence of feedback. Finally, while clouds 
can in many cases be disrupted over ~ lOMyr by shear in the 
disk, during this time the runaway collapse of the gas has turned 
an order-unity fraction of the GMC mass into stars (Fig. [T9j. As 
a result, although models without feedback can give rise Xo \ 
disks with several observed disk and GMC scaling relations, grav- 
ity alone cannot (in our calculations) regulate the structure of the 
ISM in a way that is similar to models with explicit stellar feedback. 

9.2 How Does Feedback Regulate Star Formation? 

What maintains the low star formation efficiency in the presence 
of feedback? One interpretation is that the instantaneous star for- 
mation efficiency of ~ 1 % per dynamical time can be explained as 
a consequence of the statistics of the lognormal turbulent density 
field (Nordlund & Padoan 1999: Vazquez-Sem adeni et al.|2003| [n| 
let al.,2004, .Krumholz & McKee.2005) . More recent - and much 
higher resolution - numerical work has highlighted, however, that 

This collapse might be very different in the presence of magnetic fields 
which can torque the GMCs and remove this angular momentum. 



star formation in such models proceeds more rapidly than previ- 
ously appreciated (~ 30% per dynamical time) if turbulence is 
driven solely on large scales jPadoan & Nordlund|201 1 ). In these 
calculations, the rate of star formation is indeed set by the mass of 
gas above a critical density that can collapse; however, when self- 
gravity and rapid cooling are included in simulations with sufficient 
dynamic range, self-gravitating regions can "detach" from the large 
scale flow and collapse without limit. Even if their internal instan- 
taneous SFR per dynamical time remains low, such self-gravitating 
regions quickly reach such high densities (short %„) that the global 
SFR per dynamical time becomes much larger than observed. 

This dynamical decoupling of high density regions is very 
similar to what occurs in our simulations without feedback. Recall, 
these maintain a similar large-scale velocity dispersion as models 
with feedback. [Bournaud et al.| |2010) also explicitly showed that 
high-resolution galaxy models with and without feedback maintain 
similar turbulent power spectra, but that the turbulent cascade from 
large scales is unable to dissociate clumps on small scales once 
they have become self-gravitating. This leads to large deviations 
from a lognormal density distribution (Fig. |12^ , which leads to a 
star formation efficiency that is an order of magnitude larger than 
observed. 

The difficulties of models in which turbulence is driven solely 
on large scales (via self-gravity or some other mechanism) points to 
the critical role of stellar feedback processes that can regulate star 
formation on small scales (e.g., by disrupting GMCs). Moreover, 
the high star formation efficiency found in high-resolution "turbu- 
lent box" simulations and feedback-free galaxy simulations sug- 
gests that the feedback required to dissociate dense, self-gravitating 
regions (GMCs and cores within them) regulates the SFR rather 
differently than simple homogeneous turbulence driven on large 
scales. 

Which feedback mechanisms prevent the runaway collapse 
seen in the no-feedback models and in models in which turbulence 
is driven solely on large scales? We show that this depends on the 
physical properties of the galaxy model, in particular the density 
of the ISM. In high-density regions (typical of the global ISM of 
starbursts and high-redshift galaxies and the interiors of GMCs), 
we find that radiation pressure is the dominant feedback mecha- 
nism (especially multiple-scatterings of IR photons). Lower-mass 
clouds (more prevalent in the MW and SMC models), with lower 
mean densities, can be similarly disrupted by HII photoionization 
(e.g., Matzner 2002). In Paper I we showed that radiation pressure 
alone leads to galaxy models that lie roughly along the Kennicutt- 
Schmidt relation; including the additional feedback mechanisms in 
this paper typically results in smaller (factor ~ 2 — 3) global SFRs 
by suppressing the collapse of gas into GMCs and by expelling 
more gas in galactic winds. However, these additional feedback 
mechanisms have relatively little effect on the integrated lifetimes 
or star formation efficiencies of GMCs themselves. 

In lower-density regions of the ISM (typical of MW-like sys- 
tems and dwarf galaxies), global self-regulation is much more 
strongly affected by heating processes (Fig. [7](. Removing radia- 
tion pressure in these regimes can still lead to runaway collapse of 
the dense star-forming gas. However, although dense regions run 
away, the feedback from the resulting star formation can act on the 
remaining lower-density gas to limit further collapse and help glob- 
ally regulate the star formation. Overall, SNe are the most impor- 
tant heating mechanism, since their integrated energetics dominate 
over stellar winds and photoionization-heating; however, within 
GMCs, stellar winds and photo-ionization are typically more im- 
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portant because they act on < Myr timescales while SNe do not 
explode until a few Myr later - longer than the GMC lifetime. 

9.3 Generic Consequences of Strong Stellar Feedback 

As discussed in ^9.1| we find that gaseous star-forming disks ro- 
bustly self-regulate at 2 ~ 1, which in turn sets the required ve- 
locity dispersion a and the disk scale height (Fig. |9]l. This is in- 
dependent of the details of feedback coupling mechanisms. In par- 
ticular, energy/momentum input in excess of what is needed for 
(2 ~ 1 will be "lost" (in e.g. galactic winds) or eventually dissi- 
pated (since cooling times are short), and any gas with 2^1 will 
collapse, turning into stars and triggering more efficient feedback. 

Likewise, we find that the properties of GMCs largely reflect 
the conditions for gravitational collapse starting from a.Q^ 1 disk. 
More precisely, this is true provided that feedback makes GMCs 
relatively short-lived, so that their instantaneous conditions are not 
that different from those "at formation." In this context the ~ M^^ 
mass function of GMCs is a consequence of the scale-free na- 
ture of gravity, with a cutoff near the global leans mass (Fig.|I4[l. 
The observed size-mass and linewidth-size {6v oc I^'^) relations of 
GMCs reflect conditions for collapse in a supersonically turbulent 
2 ~ 1 disk (Figs.[T5jT6j. This in turn requires a weak dependence 
of GMC densities on the mean density of the galactic disk, such 
that the linewidth-size relation normalization scales with the sur- 
face density of the disk as Sv/R°'^ oc E'' ^ 

The fact that the GMC properties - and some of the average 
disk properties - are largely determined by the conditions for grav- 
itational collapse accounts for why other calculations with different 
physics reproduce many of the same GMC/disk properties. The key 
criterion is that feedback is sufficient to both disrupt GMCs and 
maintain the ISM in a marginally stable state against collapse. If 
this can be satisfied, galaxy-averaged SFRs are set by the number 
of young stars needed to supply the necessary feedback, indepen- 
dent of how those stars form. As a result, the global SFR - and thus 
the global Schmidt-Kennicutt law - emerges relatively independent 
of the assumed micro-scale (high density) star formation law, den- 
sity threshold, or efficiency. In Paper I we varied these parameters 
widely and explicitly showed they have no effect on the SFR. Here, 
we similarly find that explicit models for molecular chemistry in 
the local SFR and cooling function have no effect on the SFR for 
the galaxies we model (Fig.[8|. 

In low-resolution numerical simulations, an explicit treatment 
of the molecular fraction can make large differences to the SFR, 
because such models are really implicit models for the fraction of 
gas at high densities. But in our calculations GMCs are resolved, 
and star formation occurs only at the highest densities within them 
(n S> lOOcm^''), which are already overwhelmingly molecular. At 
our resolution, in fact, most molecular gas is non-star forming. As 
long as a cooling channel is available, gas will eventually collapse 
to the densities needed to self-shield and form sufficient young stars 
to regulate the galaxies' SFR via feedback. An explicit treatment of 
the molecular chemistry is thus only dynamically important at low 
enough metallicities that the cooling time by any channel becomes 
long relative to the dynamical time, roughly Z ^ O.OI Zq (see also 
[Glover & Clark|20I2| ) 

These conclusions are reassuring in the sense that they suggest 
that our models are unlikely to be inaccurate on large scales, even 
if they are still missing some important physics. For example, mag- 
netic fields and/or cosmic rays may contribute comparably to the 
ISM pressure in many galaxies. However, given that our results are 
robust to variations in the "microphysics" of how, precisely, feed- 



back and cooling are implemented - and self-regulation at 2 ~ 1 
implies they will, in equipartition, contribute only ~ 10 — 30% level 
corrections to the SFR needed to maintain stability - they are un- 
likely to significantly alter our conclusions. 

The phase structure of the ISM is much more sensitive to the 
microphysics of stellar feedback than global disk or GMC prop- 
erties. Without feedback, there is a runaway collapse of some of 
the gas at the highest resolvable densities n ^ 10''cm~\ even 
though the amount of gas in GMCs and at modest densities n ^ 
1 — lOOcm"^ is similar to models with feedback (Fig.|l2|. Based 
on the calculations in this paper and in Paper I we conclude that the 
fraction of the ISM mass at high densities varies continuously with 
feedback strength, with more efficient feedback implying less gas 
at high densities (e.g.. Fig. 8 of Paper I). This is because the SFR it- 
self adjusts depending on the efficiency of feedback and the SFR is 
set by the amount of mass in the high density ISM. The dense gas 
is particularly sensitive to radiation pressure: removing this sup- 
port allows GMCs to collapse to much smaller radii and higher 
densities, particularly in high surface density systems where pho- 
toionization feedback is less effective (Fig.[6]&|15^. These results 
strongly suggest that observational probes of high density molecu- 
lar gas (e.g., HCN) are particularly sensitive probes of the dominant 
stellar feedback mechanisms in galaxies. We will pursue this con- 
nection quantitatively in future work. 
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APPENDIX A: NUMERICAL TESTS OF THE 
LONG-RANGE RADIATION-PRESSURE MODEL 



0.04 



The model we introduce in § |2.1.5| for the treatment of long-range 
radiation pressure forces requires several simplifying assumptions, 
in lieu of the ability to perform full on-the-fly radiation hydrody- 
namics with dust. In this Appendix, we test these assumptions in 
detail to calibrate our simplified model. 



Al Calibration of tlie Model 

First, we wish to check whether our calculation returns a reason- 
able approximation to the correct column densities and extinction 
of each star particle, and the corresponding rate at which the stellar 
momentum couples to the gas. 

We estimate the line-of-sight extinction of sources with our 
local column density estimator defined in Equation [6] which uses 
the local density and density gradient. If the column is dominated 
by gas close to the sources, this should be a good approximation to 
the true column towards incident gas particles. Figure [ATj compares 
this estimator to the exact line-of-sight column density integrated 
between each source and every gas particle in the simulation. We 
calculate the latter in post-processing at a given instant, integrating 
over the full simulation gas distribution. The top panel in Figure [AT| 
compares the distribution of columns towards all sources, while 
the bottom panel shows the fraction of column density ratios (local 
estimate/true) averaged over all sources. The results Figure [AT] are 
for our "standard" (all feedback enabled) HiZ model at a random 
time once the disk is in quasi steady-state, but the results are similar 
at all times and in the other simulations (modulo the absolute value 
of the typical columns). We focus on the HiZ case because it is the 
simulation where the radiation pressure has the largest effect, and 
therefore is most sensitive to any errors in the model. 

The two methods yield quite similar column density distribu- 
tions. The local estimator yields a slightly broader and more uni- 
form distribution, which is expected since it introduces some er- 
rors that broaden the distribution. This is evident in the distribution 
of A'// (estimated) /A'h (true), which has a relatively narrow Gaus- 
sian core. The dispersion in that distribution is reassuringly small 
- a factor of ~ 2 - and most important the mean value is about 
unity. The biggest difference appears at very low columns: below 
Nh ~ 10^°- 10^' cm , there is a comparable column contributed 
from passing through almost any part of a typical galaxy disk, so 
the "true" columns are likely to have contributions from large radii 
that are not accounted for in our local estimator. This is manifest 
as the non-Gaussian "tail" towards lower A'H(estimated)/A'H (true). 
However only a relatively small fraction of the stellar luminosity 
(~ 10 — 20%) samples this tail in the distribution. 

If we consider all sources by age separately, the youngest stars 
tend to have the highest local column densities, and therefore be 
most accurately approximated by our local column estimator. This 
is because they tend to be located in the sites of star formation 
(near dense gas) and often have not cleared away all of the sur- 
rounding gas and dust with feedback. If we repeat the comparison 
in Figure 



Al 



for "old" stars with age > lO' yr (> 10** yr), we find 
a median (estimated) /Wh (true) of ~ 0.7 (~ 0.5), a scatter of a 
factor ~ 4 10). This is because the older stars, having cleared 
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Figure Al. Top: Distribution of column densities between all sources (stars) 
and absorbing gas particles in a standard (all feedback enabled) HiZ sim- 
ulation. We compare the columns estimated using the local on-the-fly esti- 
mator used in real time in the code (based on the local density and density 
gradient) and that from a full (post-processing) integration over the gas dis- 
tribution between all stars and gas. The differences in time are minimal. The 
HiZ model is shown because it is most sensitive to radiation pressure, but 
the other models are qualitatively similar Bottom: Distribution of the ratio 
A'H(estimated)/A'H(true) averaged around all sources (weighted by their 
contribution to the bolometric luminosity). The local estimator we employ 
provides a good approximation to the true columns, with an accurate mean 
and factor ~ 2 scatter This is because columns tend to be dominated by gas 
near each star The largest errors occur at the lowest columns, which tend 
to have comparable contributions from gas at much larger radii (this gives 
the "tail" in the lower panel), but these tend to be older stars that are no 
longer in gas-rich star-forming regions, and so contribute only ~ 10% of 
the luminosity. 



away the dense gas from which they formed, can have larger con- 
tributions to the column from gas at very large distances. In prin- 
ciple, we could adopt an age-dependent "fudge factor" to account 
for this, but this introduces a number of additional uncertainties. In 
practice, since the older stars contribute fractionally less to the total 
luminosity especially in the UV, which is the band most sensitive 
to the column density estimator, this does not introduce much bias 
in the luminosity-weighted Nh distribution. As a result, the effect 
on the actual momentum coupled is small. 

Figure [A2| show s the consequences of our local estimator for 
the emergent SED and the typical SED seen by each gas parti- 
cle. First, consider the emergent spectrum. Using the SUNRISE 
code (lonsson 2006), we can calculate the emergent spectrum from 
the galaxy along a large number of sightlines, uniformly sampling 
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Figure A2. Top: Comparison of the emergent SED estimated with the on- 
the-fly simulation approximations and a full radiative transfer calculation 
using the SUNRISE code. For the same HiZ snapshot as Figure |A1| we 
plot the full emergent spectrum (plotted as \Lx) calculated with RT (solid 
lines; we show 10 of random sub-sampled sightlines). For comparison, the 
spectrum without any attenuation is shown as the black dotted line. Inte- 
grating over our three wavelength intervals (UV/Optical/IR), the resulting 
integrated luminosity is shown (open black points; circle/square/triangle for 
UV/Optical/IR, respectively) with horizontal error bars showing the wave- 
length range of integration and vertical error bars showing the it 1 cr range 
across all sightlines. Red filled points show the average luminosity in each 
band estimated using the simulation approximations. Shaded blue areas 
show the it 1 (T observed range of luminosities in each band, for observed 
galaxies with the same bolometric luminosity. The local approximation pre- 
dicts fluxes in each band within the 1 a range of sightlines from a full RT 
calculation; both agree well with the observed "average" SEDs of galax- 
ies with similar SFR. Bottom: Distribution of coupled momentum (specifi- 
cally, distribution across all gas particles of the radiation pressure momen- 
tum absorbed relative to the incident flux, which depends on the SED shape 
and opacities). We compare the simulation approximation and results again 
from a full RT calculation. The two distributions agree to within ~ 10%. 

the sightlines to e.g. gas particles outside the disk. We do so us- 
ing the exact same assumptions as in our simplified model about 
the scaling of dust-to-gas ratio with metallicity and emergent SED 
from all star particles as a function of age and metallicity from 
STARBURST99. Figure |A2] shows these results (as well as the un- 
processed spectrum) for a random subset of ten sightlines. Again, 
this is prohibitive to calculate in real-time, so we do so in post- 
processing for the same snapshot considered in Figure [AT| Next, 
recall that we simplify our calculation in the radiation pressure 
by integrating over three frequency ranges Luv (A < 3600 A), Lopt 
(3600 A< A < 3/.t) and Lir (A > 3^). We therefore integrate the 
emergent SUNRISE spectra over each range, and plot the result- 



ing median luminosity in each range as well as its dispersion (here 
14-86% intervals). 

Figure [A2] also compares this 'exact' SED calculation to the 
approximate result given by our local estimator. Specifically, using 
the results from STARBURST99 to estimate the initial Luv /iopt 
for each star particle, we attenuate each according to the local col- 
umn density estimator calculated for the particle (eqn. |5j and as- 
sume all of the absorbed energy in each of the UV and optical 
ranges is re-radiated as Lir. The resulting SED agrees well with 
the full radiation transfer calculation; at all wavelengths it is within 
the 1 (7 sightline-to-sightline scatter of the full calculation, and gen- 
erally within a factor of < 2 of the mean in each waveband. In the 
IR, which contains most of the energy, the accuracy is much better, 
^ 10% (the largest uncertainties come when the net attenuation is 
very large). It is important to compare these theoretical fluxes to 
those measured in similar observed galaxies. In § |A3| below, we 
discuss a variety of observational constraints on the fraction of the 
luminosity observed to escape in the UV/optical/IR, for galaxies of 
similar bolometric luminosity and/or SFR to the HiZ model consid- 
ered here. From the references discussed there, we estimate an av- 
erage empirical SED shape (shaded blue bands in Fig. |A2| (, which 
is reasonably consistent with the model predictions. 

Given the agreement in the emergent SED, we would expect 
the local model to reasonably capture the stellar momentum that 
couples to the gas. Figure |A2| tests this explicitly. For the same 
simulation, we use the full radiative transfer results to calculate (in 
post-processing at a given instant) the net momentum deposition 
from radiation absorbed in each gas particle (using the full spec- 
trum and dust opacity curve). We compare this to the actual cou- 
pled momentum using the simulation code approximations. Since 
we already know the total photon momentum is the same in either 
case, is is convenient to normalize the absorbed momentum by the 
incident momentum flux. This then tests whether the differences in 
emergent spectral shape (from our local column estimator), and the 
approximation of three wavebands instead of a full spectrum, as 
well as neglect of re-radiation changing the spatial distribution of 
emission on kpc scales, ultimately make a difference to the coupled 
momentum. Figure [A2] shows that the differences are small, at the 
level of tens of percent. This is smaller than the systematic uncer- 
tainties in our method (e.g., dust properties, unresolved clumping). 

A2 The Effects of Photon "Leakage" 

One subtle complication is that the ISM is probably inhomoge- 
neous on small scales well below what we model. In principle, that 
can allow photons to "leak out" from the region around their parent 
star at a rate much higher than the nominal exp(— ro) expectation 
which we use above to estimate the attenuated spectrum. In Paper 
I we discuss how this leakage in the IR can affect the short-range 
radiation pressure term (and show the effects are generally small 
for the cases modeled here). Here, the question is how this might 
modify the initial attenuation of optical and UV photons used to 
treat the long-range radiation pressure forces. 

If a medium with mean opacity ro has a "true" column density 
distribution /'(T|ro) across all sightlines, then the escape fraction 
is just 

/e.c = I exp (-r) dPij 1 ro) . (AI) 

Fortunately, in Paper I we show that we can make a reasonable es- 
timate of dP{T I ro). In ultra-high resolution simulations, we calcu- 
late dP{T \ ro) with ~ 1000 lines-of-sight for each molecular cloud. 
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Figure A3. Effects of allowing for photon "leakage" - i.e. assuming that 
the ISM is clumpy and porous on sub-resolution scales. We plot the ratio of 
coupled momentum from the UV and optical wavebands - the distribution 
of this quantity is measured over all gas particles for a given instant in the 
same HiZ simulation. The "standard" model takes the local column density 
estimator as-is, and estimates the escape fraction from the spectrum of each 
star with /esc = exp (— r) for each waveband. The "with leakage" model as- 
sumes that each star is actually surrounded by an (un-resolved) distribution 
of columns following Equations |AHA3| and uses this to calculate /esc. The 
allowance for leakage only makes a difference when the average opacity 
is large, at which point the contribution (in either case) to the global total 
luminosity in the relevant band is small. As a result it has a negligible effect 
on our results. 



and find the resulting column density distribution (per cloud) can 
be well-approximated with a near-universal function 



dP{r\To) 
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with a = 0.25 — 1.0 (median 0.5). This is also similar to the distri- 
bution estimated in even higher resolution simulations of individual 
GMCs and sub-cloud clumps (see'Ostriker et al."2001\ and to ob- 
servational estimates jWong e t al. 2008; Goodman et al. 20091. If 
the ISM is self-similar, then for any given to the true r distribution 
is given by the function above; we then solve for /esc: 



,^esc = ^ [vE, + ,/,(ro) + r(i)-r,/,(ro)] 



(A3) 



where the E and F are incomplete exponential integral and Gamma 
functions. The exact functional form, however, is not important: the 
behavior is simple in two limits. For ro <C 1, this function is very 
close to exp(— To), since most sightlines are optically thin (actu- 
ally /esc — > exp (—To /[I — a"]), a slightly lower value, owing to the 
presence of some dense sightlines). For ro 3> 1, 



ir[l + a-']ro-'/" 



(A4) 



i.e. the escape fraction declines as a power-law oc Tq"''''^, rather 
than as an exponential exp (—to). This is straightforward to under- 
stand - /esc in this limit is essentially just the fraction of sightlines 
which are optically thin. 

Figure A3 shows the effects of replacing /esc ~ exp (—to), our 
usual estimator of attenuation, with this modified /esc, in the HiZ 
simulation. We specifically show the distribution of coupled mo- 
mentum on all gas particles from the UV and optical wavebands 
There is no significant difference. The reason is simple. For inter- 
mediate or low To, the escape fractions are large and very simi- 
lar regardless of sub-structure. For high to, they can be very dif- 



ferent, but they are both still very small. At to ~ 100 (common 
for the optical opacities in dense star-forming regions), for exam- 
ple, the power-law /esc is much larger than exp (—to), but is still 
~ 10^' — lO^'' (for a = 0.25 — 0.5), so the amount of escaping op- 
tical light is still totally negligible compared to gravity. In short, the 
total UV/optical light seen by a given gas particle tends to be dom- 
inated by the stars that have low to (even if they are rare), not by 
the small leakage from stars with high to. For the same reasons, re- 
peating the experiment above but with a log-normal distribution as- 
sumed for dP{T I To) (instead of the power-law distribution above) 
also makes no difference. 



A3 Comparison to Purely Empirical SED Models 

The average SED shape as a function of galaxy star formation rate 
or bolometric luminosity is well-studied (see e.g. [Cortese et al. 



2006; Iglesias-Paramo et al. 2006 i|Buat et al.||20 07; Redd y et al. 
12008; Sarg syan & Weedman|200^|Noll et al.'l2009||Ta keuchi et ak 
pOlO). Specifically, using the relations inlCortese et al. (20061 be- 
tween Luv, LiR, and //-band luminosity, with the same templates 
discussed therein to make the (small) correction between their 
wavebands and ours, we expect a typical luminosity fraction in 
each interval (/uv, /opt, /ir) equal to (0.05, 0.2, 0.75) for the HiZ, 
(0.13,0.40,0.47) for the MW, (0.3,0.3,0.4) for the SMC, and 
(0.07,0.23,0.7) for the She models. We stress that all of these val- 
ues have factor > 2 scatter even in narrowly-selected galaxy sam- 
ples. This SED is shown for the HiZ model in Figure [A2l Although 
the C ortese et al.| ( |2006^ sample is low-redshift, a similar result is 
obtained if we consider the high-redshift relation studied in [Reddy] 
et al. ( 2008 1. The shaded range in Figure|A2]corresponds to the em- 
pirical scatter in the luminosity at each band at fixed L\,o\ (estimated 
as the 14 — 86% interval from the observed samples). The observed 
SED shape agrees well with the integrated fractions predicted by 
our models and the more detailed radiative transfer calculation. We 
find similar agreement for each of the other galaxy models as well. 

One way we can test our long-range radiation pressure model 
is to calculate the bolometric luminosity and incident force on gas 
as in § |2.1.5| but to specify the SED shape of each star particle 
(the emergent L\jy /Lopt/Lm), to be a constant determined by the 
observed integral values for each galaxy model. This will of course 
gloss over spatial and temporal differences in the SED shape, but it 
has the advantage that we can directly take the SED motivated by 
observations for all "typical" particles. 

In Figure [A4| we compare the star formation history for each 
of our standard galaxy models to such a calculation with a fixed 
attenuated SED. It is reassuring that the results are similar. Of 
course, this is not that surprising given that the predicted SED 
(.Ajv, ,/opt, /ir) from our simple attenuation model is similar to 
those observed in each galaxy type. Moreover, since the acceler- 
ation by the diffuse stellar radiation is primarily important only for 
gas that is out of the midplane, it turns out not to be a very bad ap- 
proximation to assume that the typical particle sees much the same 
spectrum that the typical observer would see. 

For contrast, we also compare in Figure [A4] the results if we 
were to enforce a uniform SED shape with a higher UV/optical 
fraction than is actually observed - which will substantially boost 
the coupled momentum. For the HiZ, Sbc, and MW models, this 
serves to dramatically suppress the SFR (especially at late times, 
because more of the gas is blown out of the disk in a wind). The 
SMC however is sufficiently low-luminosity that this has only a 
weak effect. This comparison highlights at least two important 
points. First, on large scales, systems are generally not optically 
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Figure A4. Comparison of our standard long-range radiation pressure model (where the attenuation and SED shape of each star particle is calculated) 
to one in which we simply force the SED shape to match a specific empirical template. We show the resulting SFR versus time (as Figure |7j for each 
of the HiZ/Sbc/MW/SMC models. In each, we compare our standard model (black) in which the SED shapes are self-consistently calculated on-the- 
fly, to a model where the SED shape (Lboi is still determined self-consistently) is forced to a constant fixed value chosen to match the observed mean 
SED shape for similar observed galaxies to each model (blue). These choices coirespond to a relative proportion of Lbol in the UV/optical/lR bands of 
(0.05, 0.2, 0.75), (0.07, 0.23, 0.7), (0.13, 0.4, 0.47), (0.3, 0.3, 0.4) for the HiZ/Sbc/MW/SMC cases. In each case, the self-consistent and empirically fixed 
models give SFRs in reasonable agreement; this is because the emergent SEDs calculated with the full models tend to agree well with the typical observed 
values (Figure [A2) . However factor < 2 differences do arise, largely because the empirical model does not allow spatial and/or time-dependent variations in the 
emergent SED (for example, very young stars in the galaxy nucleus tend to be much more obscured than older stars in the galaxy outskirts). We also compare a 
model with fixed UV/optical/IR fractions, but much higher UV/optical fractions than actually observed (red): (0.25, 0.3, 0.45) for the HiZ/Sbc/MW cases and 
(0.4, 0.4, 0.2) for the SMC. This artificially boosts the long-range radiation pressure, and significantly suppresses the SFR in the HiZ/Sbc/MW models (the 
effect is weak in the SMC case, because radiation pressure is relatively unimportant). Accounting for attenuation changing the SED shapes and consequently 
radiation pressure is critical; but small differences between empirical values and our simulation values are not dominant uncertainties in the model. 



thick in the IR and so, unlike inside of GMCs, it is the UV/optical 
photons that tend to be absorbed and provide the radiation pres- 
sure force, so it is important that any model be able to reason- 
ably approximate the average attenuation around stars in a self- 
consistent manner. Second, if there are times when dust obscura- 
tion is relatively weak, for example at very high redshifts (owing 
to low metallicities) or perhaps after phases of strong AGN feed- 
back (either the AGN blowing away or dissociating the dust), then 
the long-range radiation pressure from stars can be dramatically 
boosted (even while the short-range effects are weakened), enhanc- 
ing global outflows. Since this in turn removes more material, there 
is the possibility of unstable runaway that could substantially en- 
hance quenching in such systems. 



APPENDIX B: RESOLUTION & CONVERGENCE TESTS 

In Figure [BT] we perform a spatial and mass resolution test for our 
HiZ simulation. We compare low, intermediate, standard, and high- 
resolution simulations of otherwise identical HiZ runs with all feed- 
back enabledf^We also compare, at low resolution, two additional 
runs with the same mass resolution but force resolution varied by 
a factor of ~ 5. The star formation history converges reasonably 
quickly; the differences at the peak of the SFR are small even at 
low resolution. This is expected and discussed in detail in Paper I; 
since the equilibrium SFR is feedback-regulated, set by the balance 
of momentum injection with dissipation, it depends primarily on in- 
tegral quantities (the total momentum injection oc L oc M, , approxi- 
mately). As a result, it converges relatively quickly. In all cases, the 
higher-resolution runs begin to rise in SFR (and generally converge 
to their equilibrium structure) slightly earlier - this simply owes to 
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Figure Bl. Resolution tests. We plot the SFR vs time as Figure|7] We con- 
sider a "standard" (all feedback enabled) HiZ model at a series in mass and 
force resolution (Wgas is the number of gas particles in the star-forming disk, 
which determines the mass resolution; egas is the force softening). The SFR 
converges quickly because it depends on the integral constraint of momen- 
tum feedback (oc L oc A/* ) balancing gravity (differences at late times owes 
to more efficient winds at higher resolution). Recall, the Jeans length in 
these galaxies is several hundred pc (Jeans mass > 10* Mq), so all of these 
models formally resolve the Jeans scales. At lower resolution, our presciip- 
tions do not have clear physical meaning. 

the fact that the most rapidly collapsing spatial scale resolved is on 
a smaller timescale as we go to higher resolution. At late times, the 
higher-resolution simulations have a slightly lower SFR; this owes 
to their blowing slightly more efficient galaxy-scale outflows. 



For technical reasons to ensure quantities such as random number gen- 
eration were identical in these runs, this set was run on the same set of pro- 
cessors and so does not overlap with our standard HiZ run shown earlier; 
but the results from that run are completely consistent with the standard- 
resolution case here. 



© 0000 RAS, MNRAS 000, 000-000 



